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Oscillating Triangles 








These figures show one of the several remarkably stable and robust structures of the equations that model double diffusion 
convection. In such flows, thermal effects and density stratification (salinity) compete to form complex equilibrium structures. 
The snapshots in this figure show the cross section of the state of a periodic solution known as “oscillating triangles " where 
the structure is neither a standing wave nor a traveling one. For more details of the possible patterns in the double diffusion 


equations, see the article by Y. Renardy in this issue 
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Dynamical Systems and 
Ocean Dynamics 


Reza Malek-Madani, Guest Editor 
Office of Naval Research 


Introduction 


Since the seminal works of Liapunov and Poincaré, re- 
search in dynamical systems has provided scientists with a set 
of practical tools to investigate the qualitative behavior of 
solutions of differential equations. In the past fifty years this 
area of mathematics has produced some of the most insightful 
and precise information about the nature of the dynamics 
present in many important applications ranging from biology 
to physics. Recently, a group of applied mathematicians 
joined together in a concerted effort to develop dynamical 
systems tools appropriate for problems in ocean dynamics. 
The articles that follow will give an introduction to these 
methods and provide a sample of the type of questions that one 
can properly pose and answer in a dynamical system setting, 
and point to the future directions of research. 

The mathematical modeling of the ocean is a relatively 
young science that still takes many of its cues from observation 
and experimentation. Although the basic governing equations 
have been well-known for many decades, a large class of 
partial differential equations are present in the literature whose 
derivations depend on the special scales in the geometry of the 
region of interest as well as the strength of the forces involved. 
In addition to the information one draws from these models, a 
substantial amount of data is collected, either remotely or 
in-situ, that has contributed significantly to our understanding 
of the basic phenomena as well as challenged some of the 
underlying assumptions that have led to the mathematical 
models. One approach to reducing the gap between the sophis- 
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tication and physical realism of the mathematical models and 
the information gathered in real data is to develop mathemati- 
cal tools that are not so heavily model dependent that they lose 
their significance when one mathematical model is discarded 
in favor of another. With this point in mind, the methods and 
techniques outlined in the next three articles in this journal aim 
at understanding the basic issues that describe transport and 
mixing in the ocean and atmosphere as well as quantifying the 
competition among the forces that lead to double diffusion 
convection. 

Computational models have been traditionally the widely 
used tool for predicting the response of the ocean and the 
atmosphere to various mechanical and thermodynamical proc- 
esses. In many respects, these efforts have closely followed 
the development of the numerical analysis of the Navier- 
Stokes equation and have attempted to find the appropriate 
range of various viscosity parameters that would lead to rea- 
sonable numerical simulations. This strategy has paid off 
well in the context of low Reynolds number dynamics mod- 
eled by this equation while it is a nich area of research for high 
Reynolds number turbulence. Because the Navier-Stokes 
equations are generally considered the correct mathematical 
model for many applications, it is felt that the investment of 
resources in developing more accurate numerical schemes for 
this specific equation is worthwhile. (See reference | for a 
viable alternative to a direct discretization of the Navier- 
Stokes equations, whether using finite difference methods, 
finite elements, or the various quasi-spectral methods that are 
currently being used in the literature). It is not clear that this 





strategy is correct, or even desirable in ocean dynamics, con- 
sidering the fact that more than one candidate for the “stand- 
ard” model exists. It then seems reasonable that the solutions 
of the governing differential equations are understood at a 
level beyond the usual listing and tabulating the various values 
of the parameters and instead search for the coherent structures 
that are commonly present in these models. 

This is not to diminish the role of numerical analysis in 
ocean dynamics. In fact, quite the opposite result is antici- 
pated. The intent of the new dynamical systems methodology 
is to elucidate the common complex structures in the set of 
relevant systems of partial differential equations. Once the 
coarse properties of these structures are identified, numerical 
simulations on a reasonably small and computationally tracta- 
ble number of initial conditions will reveal the finer properties 
of the mathematical models. 


Scope of Applications 


The mathematics described in the following articles is 
divided in two categories, depending on whether the models 
have a Lagrangian formulation or an Eulerian one. In the 
articles of S. Wiggins” and C. Jones et. al.’ of this issue, the 


goal is to understand the nature of the solutions of the per- 
turbed systems of ordinary differential equations that arise in 
the Lagrangian setting of motions in the ocean and atmos- 
phere. This point of view is specially significant in ocean 
dynamics in light of the type of data taken from the so-called 
Lagrangian drifters that provide real data about various fluid 
particles in the ocean. Because only a finite (and a relatively 
small) number of these drifters are available at the present 
time, one is left with the challenging task of gathering infor- 
mation about a portion of the ocean from a handful of particle 
trajectories. This fact alone shows the acute need for a kind 
of mathematical analysis that goes well beyond the tabulation 
of parameter values. 

The methods that are discussed in references 2 and 3 show 
how certain structures associated with invariant regions of an 
unperturbed dynamical system are transported by the per- 
turbed velocity fields. Specific and delicate mathematical 
tools are developed to monitor the transport of information for 
a rather large class of perturbations. One of the intriguing 
problems that will be addressed in the future by these investi- 
gators is what happens to these dynamical systems methods 
when, as is the case in real life, the velocity fields that the 
differential equations are based on are only known on a dis- 
crete set of points in the ocean. A related problem is whether 
one can develop systematic mathematical tools that can dis- 
tinguish the effect of chaos from stochasticity in data sets that 
exhibit complex and coherent structures. 

The article by Y. Renardy of this issue* addresses the set 
of partial differential equations that model double diffusion 
convection. This paper shows how to use the center manifold 


theorem of dynamical systems and classify the stability prop- 
erties of certain important solutions of the original system. A 
central feature of this theorem is in its ability to reduce the 
analysis of the dynamics to a lower order system. Various 
traveling wave solutions of the original governing system of 
partial differential equations are identified and numerically 
simulated. 
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Lagrangian Transport in 
Geophysical Flows: 
New Approaches and 
Techniques from 

amical Systems Theory 


Stephen Wiggins 


Applied Mechanics, and Control and Dynamical Systems 


California Institute of Technology 


Introduction 


The standard map is a two-dimensional area-preserving 
transformation, or dynamical system, that has been the subject 
of literally hundreds of papers in the dynamical systems com- 
munity over the past twenty years. The map is given as follows 


x—> xX+y-T, 
y — y+ksin(x+y), 


and it is periodic in x and y. In Figure 1 we plot a single orbit 
of this map for different initial conditions in each of the four 
panels for the parameter value k=1. (By the term “orbit” we 
mean a sequence of points obtained by iterating a given initial 
point according to the map. Thus an orbit of a map is a discrete 
time trajectory.) The orbits in the upper left hand and lower 
right hand panels explore a significant fraction of the phase 
space. Moreover, the motion is chaotic in the sense that nearby 
orbits separate at an exponential rate. However, note the 
“holes” in the regions. The chaotic regions do not fill all of the 
space. Note the difference in the orbits in the upper right hand 
and lower left hand panels. These orbits appear to lie on closed 
curves. There are four discrete closed curves in the former case 
and one closed curve in the latter case. Orbits on these curves 
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separate at a linear rate. These closed curves are examples of 
“KAM tori” and are boundaries which orbits cannot cross. 
KAM tori are examples of lower dimensional “phase space 
structures” that evidently have a marked influence on transport 
in phase space. 

The chaotic regions are built on the skeleton of a lower 
dimensional phase space structure. In Figure 2 we show some 
intertwining one dimensional curves that lie in the chaotic 
region in the upper left hand panel of Figure 1. These curves 
are the one dimensional stable and unstable manifolds of the 
saddle type fixed point in that region. They are the mechanisnt 
giving rise to the chaos, and they also govern the transport 
through such chaotic regions and play a very important role in 
understanding this phenomena. 

Thus, in this example we see that the phase space is a 
mixture of “regular” and chaotic regions, a situation that we 
now know is typical for many nonlinear dynamical systems. 
Moreover, lower dimensional “phase space structures”-KAM 
tori, stable and unstable manifolds of periodic orbits—form the 
skeleton in phase space upon which the regions of regular and 
chaotic mouons exist. 

Despite the simplicity of this example, the same mathe- 
matical framework and ideas have implications and uses for 
many problems in geophysical fluid dynamics, and we next 





want to motivate this setting and approach. Suppose one is 
interested in the motion of a passive tracer in a fluid (e. g. dye, 
temperature, or any material that can be considered as having 
negligible effect on the flow), then, neglecting molecular 
diffusion, the passive tracer follows fluid particle trajectories 
which are solutions of 


x=v(x,t;p), (1) 


where u(x,t) is the velocity field of the fluid flow, 
x € IR", n=2 or 3,and p € JR? represent possible parameters. 
When viewed from the point of view of dynamical systems 
theory, the phase space of (1) is actually the physical space in 
which the fluid flow takes place. Evidently, “structures” in the 
phase space of (1) should have some influence on the transport 
and mixing properties of the fluid. To make this more precise, 
let us consider a fluid mechanically more simplified situation 
that can be related directly to the standard map example 
described above. Suppose that the fluid is two-dimensional, 
incompressible, and inviscid, then we know that the velocity 
field can be obtained from the derivatives of a scalar valued 
function, w(x, ,x,,;11), known as the stream function, as follows 


= 
x,= x, (X, X45), 


> (x, X61), (x, .x,) €IR?. 





Figure 1. 


Four orbits of the standard map for k=1. 




















In the context of dynamical systems theory, (2) is a time-de- 
pendent Hamiltonian vector field where the stream function 
plays the role of the Hamiltonian function. If the flow is 
time-periodic then the study of (2) is typically reduced to the 
study of a two-dimensional area preserving Poincaré map, of 
which the standard map is an example. Practically speaking, 
the reduction to a Poincaré map means that rather than viewing 
a particle trajectory as a curve in continuous time, one views 
the trajectory only at discrete intervals of time, where the 
interval of time is the period of the velocity field. The value of 
making this analogy with Hamiltonian dynamical systems lies 
in the fact that a variety of techniques in this area have 
immediate applications to, and implications for, transport and 
mixing processes in fluid mechanics. For example, the persist- 
ence of invariant curves in the Poincaré map (KAM curves) 
gives rise to barriers to transport, chaos and Smale horseshoes 
provide mechanisms for the “randomization” of fluid particle 
trajectories, an analytical technique, Melnikov’s method, al- 
lows one to estimate fluxes as well as describe the parameter 
regimes where chaotic fluid particle motions occur, a rela- 
tively new technique, lobe dynamics, enables one to efficiently 
compute transport between qualitatively different flow re- 
gimes, bifurcation theory enables one to understand how 
qualitatively different flow regimes appear and disappear as 
system parameters are varied, etc. In short, many new methods 
for the study of fluid transport and mixing are available from 
dynamical systems theory. In the following we will show how 
this framework and techniques can be used in a problem of 
geophysical interest. 

The techniques and results described in the following 
have been obtained over the past few years in collaboration 
with my colleagues Anthony Leonard and John Brady at 
Caltech, and former graduate students Vered Rom-Kedar, 
Roberto Camassa, Tasso Kaper, Darin Beigie, and Igor Mezic’. 


Transport in a Cellular Flow: 
The Physical Setting and the 
Model Flow 


The physical setting that we consider is as follows. We 
consider a (steady) convection cell whose horizontal length is 
much larger than its height and where the convection cells are 
aligned along the y-axis (Figure 3). In this situation the flow 
is essentially two-dimensional and, assuming stress-free 
boundary conditions and single-mode convection, an explicit 
form of the velocity field is given by (see Chandrasekhar 
(1961]) 


x=- Az cos(1z) sin(kx) = 2 (x,2), (3) 
k dz 


z = Asin(nz)cos(kx) = Mes ; (4) 
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wix2) =4 sin(kx) sin(n), 


A is the maximum vertical velocity in the flow, 


2n 
ial | 


(A the wavelength associated with the cell pattern), and the 
length measures have been non-dimensionalized so that the 
top is z=1 and the bottom z=0. This flow has a countable 
infinity of saddle type stagnation points, or hyperbolic fixed 
points in dynamical systems terminology, on the upper bound- 
ary at 


(x,2) = a )j =0,41,42.... 


and a countable infinity of hyperbolic fixed points on the lower 
boundary at (z,2)= En, = 0+1,#2,... Fixed points with the same x 


coordinate are connected by a stagnation streamline, or het- 
eroclinic orbit in dynamical systems terminology. The result 
is an infinite number of cells, and many of the fluid mechanical 
problems of interest are concerned with the transport of a 
passive tracer from cell-to-cell. Fluid flows which are divided 
up into a collection of cells occur in a variety of geophysical 





Figure 2. 


Stable and unstable manifolds of a period one point in a chaotic 
region. 
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applications and the methods and approaches discussed in this 
article in principle can be applied to any such cellular flow. 


The Experiments of Solomon 
and Gollub 


We motivate the analysis by discussing some experiments 
of Solomon and Gollub on this system. If the temperature 
difference between the top and bottom of the convection cell 
is increased an additional time-periodic instability occurs, 
resulting in a time-periodic velocity field (Clever and Busse 
[1974], Bolton et al. [1986]). Solomon and Gollub [1988] 
experimentally studied the transport of dye in such a flow and 
they made the following observations. 


l. There was a dramatic enhancement of the “effective 
diffusivity”. 
Molecular diffusion appeared to play no role in the 
transport. 


The flux across the cell boundaries (heteroclinic orbits) 
depended linearly on the amplitude of the oscillatory 
instability and was independent of the wavelength A. 


Therefore, the transport problem is radically different in 
the unsteady case as compared to the steady case since in the 
steady case streamlines cannot cross and therefore tracer 
crosses the cell boundaries solely as a result of molecular 
diffusion. In order to study these issues Solomon and Gollub 
introduced the following model of the even oscillatory cell 
instability: 


x= ~“Fcos(nz)(sin(kx) + ekfit)coskx)] = Mx21) 


z = Asin(mz)(cos(kx) — ekfit)sin(kx)] = M20) 


where 
y(x,z,t) = 4 sin(kx)sin( nz) + efit)cos(kx)sin(1z) 


and we will take 
ft) = cos(ax) ande ~(R-R,)'”, 


where R is the Rayleigh number and R, its critical value at 
which the time-periodic instability occurs. The pros and cons 
of this model are discussed in Solomon and Gollub [1988]. 
This model will be the starting point of our analysis and, 
motivated by the experimental observations of Solomon and 
Gollub, we will consider the following four questions. 


1. | What is the mechanism for cell-to-cell wansport? 


2. Can we quantify the spreading of a passive tracer (dye) 
initially contained in one cell? 








Figure 3. 


Streamlines for the steady, cellular flow. 





























Can we predict the number of cells invaded as a function 
of time? 


What are the effects of the addition of a small amount 
of molecular diffusion? 


The Mechanism for Cell-to-Cell 
Transport 


In the case of unsteady two-dimensional flows modem 
dynamical systems theory provides us with some new methods 
and concepts that enable us to analyze the particle kinematics. 
For example, the saddle-type stagnation points may move and 
the separatrices associated with different “moving stagnation 
points” typically will not coincide to form cell boundaries in 
the flow. This situation can lead to cell-to-cell transport. Very 
roughly speaking, dynamical systems theory provides us with 
methods for studying particle kinematics in unsteady flows 
with moving stagnation points and their associated moving 
separatrices. For the case of time-periodic variability, if the 
variability is weak then it can be shown that a saddle-type 


stagnation point persists, but it oscillates periodically in time 
(with the period of the variability). Moreover, this “oscillating 
Stagnation point” (more precisely, hyperbolic periodic orbit) 
still maintains its saddle-type stability characteristic and asso- 
ciated with it are two one-dimensional periodically varying 
curves called the stable and unstable manifolds of the hyper- 
bolic periodic orbit. Trajectories on the stable manifold ap- 
proach the hyperbolic periodic orbit asymptotically as 1 — oo 
and trajectories on the unstable manifold approach the hyper- 
bolic periodic orbit asymptotically as t + — oe. These curves 
are invariant curves in the sense that particle trajectories that 
start on the curves remain on the curves for all time. This 
property of invariance implies that other trajectories cannot 
cross these curves, hence they act as barriers to fluid transport 
in a way that we will describe more fully shortly. These curves 
are the unsteady analogs of the separatrices that commonly 
arise in steady flows. However, as opposed to the steady case, 
the stable and unstable manifolds of different oscillating sad- 
dle-type stagnation points typically do not coincide to form 
complete barriers to fluid transport. Rather, at a fixed instant 
of time they may intersect in a discrete set of points. Never- 
theless, these curves still have the important kinematical effect 
of forming the boundaries between fluid particle trajectories 
that move in different directions after their encounter with the 
oscillating saddle-type stagnation point. However, since the 
stable and unstable manifolds of different oscillating saddle- 
type stagnation points need not join up (and, hence, form a 
curve of finite length in the flow), they may have infinitely 
long length and wind throughout the flow (and themselves) in 
a complicated geometrical fashion. Keeping in mind the kine- 
matical meaning of these curves as boundaries between quali- 
tatively different fluid particle motions, they therefore act as 
partial barriers to fluid transport. As a result, fluid transport 
in two-dimensional unsteady flows may be much more varied 
and complex. Dynamical systems theory gives us a way to 
make these intuitive notions precise and quantitative, which 
we now describe for time-periodic variability. 

For time-periodic variability a variety of standard dy- 
namical systems techniques and results make it advantageous 
to consider particle kinematics via the so called Poincaré map, 
which we will denote by f. That is, rather than plotting a 
particle trajectory as a continuous curve in space one plots its 
evolution at discrete intervals of time, where the interval of 
time is equal to the period of the variability. In this setting, the 
oscillating saddle-type stagnation point is manifested as a fixed 
point of the Poincaré map and the stable and unstable mani- 
folds associated with the fixed point are similarly fixed in 
space for the Poincaré map. Consequently, the situation bears 
some similarity to the steady case, but it is important to keep 
in mind that for the Poincaré map particle trajectories are 
manifested as sequences of discrete points, rather than con- 
tinuous curves. Also, as we mentioned earlier, the stable and 
unstable manifolds of different hyperbolic fixed points may 
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Figure 4. 


The behavior of the stable and unstable manifolds on the wails of 
the convection cell under the time-penodic perturbation— the 
heteroclinic tangle, numerically calculated fore = 0.1, o = 0.6, 
A=0.1. 











L 








not coincide(or “join up’’) to make a separatrix of finite length, 
rather they may intersect in a discrete number of points and 
wind throughout the flow. We demonstrate this for (5) in 
Figure 4. The stable and unstable manifolds of the different 
hyperbolic fixed points actually have infinite length, but we 
obviously can only show a finite portion of each in the figure. 
Next we show how the geometric structure associated with the 
stable and unstable manifolds of the different hyperbolic fixed 
points influence transport. 

For the Poincaré map, segments of finite length of the 
stable and unstable manifolds of the hyperbolic fixed points 
can be used to form cell boundaries. We illustrate this in Figure 
5 where the segments of finite length begin at the hyperbolic 
fixed points and end at particular intersection points of the 
stable and unstable manifolds, which we will refer to as the 
boundary intersection points. Now we can consider the ques- 
tion of fluid exchange between the different regions defined 
by our chosen boundaries that are shown in Figure 5. 

First we list two rules that points on the stable and unstable 
manifolds of hyperbolic fixed points must obey under iteration 
by the Poincaré map. (By “iteration” we mean repeated appli- 
cation of the Poincaré map to a particle or, in other words, 
discrete time evolution of a particle trajectory.) 


1. A point that is on both the stable and unstable manifolds 
must remain on both manifolds under all positive and 
negative iterations of the Poincaré map. This is because 
these manifolds are invariant manifolds. 


Points on the stable or unstable manifolds of a fixed 
point maintain their relative order (in the sense of dis- 
tance in arc length from the fixed point) under iteration 
by the Poincaré map. That is, if we consider two distinct 
points on the stable or unstable manifold then one will 
be closer than the other to the fixed point in the sense of 
arc length along the stable or unstable manifold. If we 
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then iterate both points, the same point will always 
thereafter be closer to the fixed point. The reasons for 
this are more technical and are related to uniqueness of 
particle trajectories for a given initial condition. 


With these two rules in mind let us return to Figure 5 and 
consider the preimages under the Poincaré map of the bound- 
ary intersection points, i. e. , the backwards time evolution of 
these points for one period of the variability. By rule 1 these 
points still lie on both the stable and unstable manifolds, and 
since the manifolds are smooth curves they must wind through 
each other as shown in Figure 6. In this figure we have shown 
also precisely one additional point of intersection of the mani- 
folds between the boundary intersection point and its pre- 
image. There is a technical reason for this that we briefly must 
mention. Generally, there must be an odd number of additional 
intersection points between the boundary intersection point 
and its preimage if we have uniqueness of parcel trajectories 
for a given initial condition. Without loss of generality, we can 
assume that we have one additional intersection point since all 
other cases can be recast in this form. Hence, the segments of 
the stable and unstable manifolds between a boundary inter- 
section point and its preimage form two lobes and these two 
lobes are said to form a turnstile. The turnstiles are the mecha- 
nisms governing transport between the different regions, as we 
will now argue. Consider the image of the turnstile lobes under 
the Poincaré map f (i. e. , the evolution of these points for one 
period of the flow). Using the two rules above, as well as the 
fact that(for a continuous and continuously invertible map) the 
boundary of a set maps to the boundary of its image and the 
interior of a set maps to the interior of its image, we see that 
the image of the turnstile lobes appear as in Figure 6. Thus, 
the respective turnstile lobes have “switched sides” in relation 
to the boundary (hence the term “turnstile”). Now, using rule 
2 above, it can be argued that the only points that cross the 
boundary in one iteration of the Poincaré map are those in the 





Figure 5. 


Cell boundaries for the unsteady flow(for the associated Poin- 
caré map). 




















Figure 6. 


The turnstile lobes associated with cell-to-cell transport. 
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turnstile lobes. Therefore the flux across the boundary in one 
period is simply the area of the turnstile lobe. 

In the case where the variability is weak there is a pertur- 
bation method that will enable us to estimate the lobe area—and 
therefore the flux per period—across the “broken separatrices”. 
This is Melnikov’s method for time-periodic vector fields. The 
Melnikov function is a measure of the distance between the 
stable and unstable manifolds measured along a line perpen- 
dicular to the unperturbed separatrix at the point parametrized 
by t), where t, can be thought of as an arc length parameter 
measuring distance along the unperturbed separatrix. In 
Camassa and Wiggins [1991] the Melnikov function on the 
zero phase cross section is calculated and found to be given 
by: 


M(tg) = casech(57-)sin(«o), 


Clearly the Melnikov function has a countable infinity of 
simple zeros, that correspond to a countable infinity of inter- 
sections of the stable and unstable manifolds. 

As described above, the flux of fluid between cells is 
given by the area of the turnstile lobes. In this perturbative 
setting an approximation to the flux (lobe area) is given by the 
integral of the Melnikov function between two adjacent zeros, 
i.¢., 


WL; o(1)) = w(Lo 1(1)) = esech~ +a(e’), 


where pt (L; {)) denotes the area of the lobe L; j We remark 
that because of translational symmetries, this result holds for 
all turnstile lobes. The symmetries of this problem are dis- 
cussed in Camassa and Wiggins [1991]. This approximation 
to the flux show us that the flux depends linearly on the 
amplitude of the oscillatory instability and is independent of 
the wavelength of the cell patterns A, exactly as observed 
experimentally by Solomon and Gollub. 


The Number of Cells Invaded 
as a Function of Time: Lobe 


Dynamics 


We next provide another example of how “lower dimen- 
sional structures” in the flow, the lobes, determine a large scale 
transport property of the flow. Suppose a given cell is filled 
with a passive tracer. How many cells are “invaded” by this 
tracer after some elapsed time? The answer to this question 
can be determined by the geometry of the lobe intersections 
associated with the cell boundaries. Using lobe dynamics type 
arguments and invariance of the manifolds, as well as transla- 
tion and reflection symmetries, an upper and lower bound can 
be placed on the time thai a passive scalar initially in R, first 
invades R_;. Surprisingly, all the necessary information is 
contained in two integers which we refer to as the signatures 
of the heteroclinic tangles forming the cell boundaries, 
namely: 


¢ m= the smallest integer such that *(L, (1))rL,_,(1)#0. 

e m’ = the smallest integer such that the boundary of 
F™**(L, (1) intersects the boundary of L_, , (1) in four 
distinct points 


The signatures are illustrated in Figure 7 for m= 1 and mf = 2. 


If we denote the period of the velocity field by 
a and :;’ the time necessary for tracer initially in R, to enter 
R; we have the following general result obtained in Camassa 
and Wiggins (1991): 


(jm + 1)T < ff < [i — lm’ + (m+ 1))T, 2 2). 
Therefore, the upper and lower bounds for the first invasion 


time are completely determined by the signatures m and 7’. This 
result is independent of the boundary conditions. For instance, 





Figure 7. 


An illustration of the signatures of the heteroclinictangle asso- 
ciated with the cell boundaries. 





2 





p 
ms 
! 
HL; o(1)) 
es OL, ” 




















Pz 
. PLo(1) 





One 1995 9 





numerical computations show that in the case A=0. 1, e = 0. 1 
and @ =0.6 we have wandw=3 indicating that one cell is 
invaded every three periods, for A=0.1, € = 0.1 and @ = 0.24 
we have mand =1 indicating that one cell is invaded every 
period, and for A=0. 1, ¢ =0.01 and @=0. 6 we have m and 77. = 
4 indicating that one cell is invaded every four periods. We 
remark that the signatures are very easy to compute numeri- 
cally as one needs to compute only a(usually short) finite 
length of the manifolds. Moreover, the necessary manifold 
interections are typically robust with respect to numerical 
computations. 


Quantifying the Spreading of a 
Passive Scalar: Lobe Dynamics 


The dynamics of the turnstile lobes, lobe dynamics, can 
be used to rigorously determine the amount of tracer originally 
in one cell that is contained in any other cell at a later time. 
Suppose the cell R, is uniformly filled with tracer (species S,) 
at t=0 how much of species S, is contained in R; at any t > 0? 
If we denote by T,. ; () the total amount of species S, inR; 
immediately after the n-th iterate and by a, (n) =T, (n)-T,, 
j(n-1) the flux of species S, into R; on the n-th iterate , then it 
can be shown that this last formula can be reduced to an 
expression containing areas of intersection sets involving im- 
ages of only one of the turnstile lobes L, (1) associated with 
the boundary of R,. The resulting expression is given by 


ay j(n) = T, (a) —T, j(n — 1) = (6 2 + 5; OL, o(1)) 


a-l 
+> (a( ly AD", ft) - uth pa HOY Ly (1D) 
el 


PL ja DVL oD) -2n(L1 VF" of) 
#1H(L joa OIL of) HL joa DVL 1D) (6) 
“HL je VL oD) +4(ly md", ol) 
#H(Ljo2_-pe1 FL 01D) ) {L241 Lyf) ) 


H(Lp2 DOSE ol) +4{L +2(DF "Li 0(1)))) 


where, recall, 1 (L, j(1)) denotes the area of the lobe L, i). 
Consequently, a relatively straightforward numerical compu- 
tation of a, (n) is possible. To compute a, (n) one merely 
locates the lobes indicated in (6), iterates L, o(1) (i.e. , a grid 
of points covering L, 9(1)), determines the area of intersection 
of the appropriate iterates of L, o(1) with the other lobes in (6), 
and adds up the result according to the recipe provided by (6). 
Numerically this method is extremely efficient and affords 
tremendous savings in CPU time over standard “brute force” 
methods; details are discussed in Camassa and Wiggins 
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[1991]. However, practically, this computation can only be 
performed for a finite amount of ume. In the next section we 
discuss an infinite time result. 


Dispersion 


Dispersion is a quantity that is frequently of interest in 
geophysical settings. For example, Pasmanter [1991] and Rid- 
derinkhof and Zimmerman [1992] have recently studied dis- 
persion in models of shallow tidal flows and found that it 
behaves asymptotically like . Weiss and Knobloch [1989] 
have performed a numerical study of dispersion in modulated 
traveling waves in binary fluid convection and found a disper- 
sion exponent of 1.93. As their results were numerical, they 
were only able to compute for a finite length of time. Using a 





Figure 8. 


Comparison of lobe dynamic transport results (solid) with trans- 
porttaking into effect molecular diffusivity (dashed) for the three 
cases. a) A= 0.1, m = 0.6,€ = 0.1, b)A= 0.1, w = 0.6, € = 0.01. 
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simple application of Birkhoff’s ergodic theorem we can de- 
rive a rigorous result that applies in each of these settings and 
gives conditions under which the dispersion should behave 
asymptotically (in time) like :7. As such, it may also prove to 
be a useful guide for numerical investigations. Our discussion 
will be in the context of the time-periodic cellular flow, how- 
ever the results can be applied more generally(See Mezic’ and 
Wiggins [1994b}, [1994c}). 

Referring to the velocity field (5), the displacement of a 
particle trajectory in the x direction is given by 


x(0)— (0) = J (AB coscaa(oytsin(iae) + 


ek cos(wt)cos(kx(0))} )at, 


where (x(t), z(t)) denotes a trajectory. For the application of 
Birkhoff’s ergodic theorem we need to introduce a standard 
dynamical systems “trick” that is applied to time-periodic 
vector fields. That is, we introduce the phase, call it 8, of the 
time-periodic part as a new dependent variable and append it 
to the vector field. The new component is therefore 6 = # and 
the vector field is then viewed as a three-dimensional time-in- 
dependent dynamical system. If D denotes the flow domain 
(periodic in x, bounded in z), then the phase space of this 
“enlarged” dynamical system is denoted D x S', where S' 
denotes the phase of the time periodic part. The “flow” gener- 
ated by this dynamical system is then denoted by ® (x,z,9). i. 
e. , this is the trajectory that starts at the point (x, z, 8). If we 
denote the x component of the velocity of (5) by u, then (7) 
can be succinctly written as 


x(t) — x(0) = J eo x2.0))de, (8) 


The mean square displacement or dispersion of the x 


component of (5) of an ensemble of points under the flow is 
given by 


((x(t) - x(0) - (x(t) - x(0)))’) = D, (0), 


where the average indicated by the angle brackets is defined 
as 


(()-2) = J cat) - x00) pa, 


P = p(x, z, 8) is the initial distribution of points (assumed to be 
bounded and integrable on D x S'), and ds denotes the measure 
or “volume element” on D x S': Incompressibility of the flow 
implies that the flow is “measure preserving”, which is impor- 
tant for the application of Birkhoff’s ergodic theorem. 


We are interested in determining the asymptotic behavior 
of the dispersion. We have the following calculations: 


a 
im 
bow 


[->00 


=lim ¢ (; u (@,(x,2,8)) dt — ¢ 7 u(P,(x,2,8))at )’, 


ss ed _ Af 2 
=( (tim : H (occa.0)ae— (ti, | u9,(x.2.0)\dh) ) ), 


= ( (u"(x,2,0) — ... u"(x,2,8)))?) =a. 


The passage from the first to the second line is justified 
by the fact that the function u is bounded and integrable 
on D x S!. In the second line the limit 


In the second line the limit 
time [U((2,0))dr = u'(x.2.8) 
-~ 0 


exists for all points in D x S' by Birkhoff’s ergodic 
theorem (see Arnold and Avez [1968]), with the possible 
exception of a set of \t-measure zero. This limit is the 
time average of the function w along the fluid particle 
trajectory that starts at the point (x, z, 6). Moreover, 
Birkhoff’s ergodic theorem also guarantees that this 
limit is integrable. This, together with the boundedness 
of u, implies that the quantity a defined above is finite, 
and if it is nonzero, we can conclude that the dispersion 
of the ensemble of particles in the z direction behaves 
asymptotically like ¢?. 





Figure 9. 


Contours corresponding to points having the same timeaver- 
age along trajectories for the x-component of velocity, u 
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Figure 10. 


The distribution of the average x-component of velocity. 
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The nature of the coefficient a gives some insight into the 
dynamical mechanism giving rise to t dispersion. It is easy to 
see that since the expression inside the angle brackets defining 
a is nonnegative, a=0 if and only if u* = <u*> on the support 
of p, i. e., on the set of points for which p(x, z, 8) is nonzero, 
with the possible exclusion of sets of measure zero. However, 
<u*> is a constant. Therefore, we can make the following 
conclusion. Let C < D x S' denote the support of p(x, z, 9). 
Then, if u*(x, z, 8) is not constant almost everywhere on C, 
D, (t) ? as t > 0, Now assume that the flow is ergodic. Then 
u*= <u*> = <u> almost everywhere. Therefore, a=0. So, a 
necessary condition for t? dispersion is non-ergodicity of the 
flow. 

We end with some final observations. 


Note that our result is independent of the Reynolds 
number. We are dealing solely with kinematical consid- 
erations. 


There has been much work done in the last 10 years on 
chaotic fluid particle dynamics in two dimensional, 
time-periodic velocity fields. In such situations one 
typically sees a mixture of regular and chaotic regions 
of fluid particle motions in the flow. Our dispersion 
result implies that if one takes an initial distribution of 
points that is, roughly speaking, not entirely contained 
in a Single regular or chaotic region, then the asymptotic 
behavior of the dispersion will go like t?. 


The case of a diffusive tracer has recently been treated 
by Mezic’ et al. [1994]. This work shows that the 2 
dispersion result in the non-diffusive case is important 
for studying the diffusive case in the high Peclet number 
limit. (The Peclet number is a dimensionless number 
that is typically of the form r=“, where U is a charac- 
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teristic velocity, | is a characteristic length, and D is the 
diffusivity.) In particular, for the time-periodic cellular 
flow it is possible to show that t? dispersion of the 
convective part of the problem implies a Pe” dependence 
of the effective diffusivity coefficient in the diffusive 
case. In fact, the dependence of the effective diffusivity 
on the Peclet number can be determined solely on the 
knowledge that the convective part of the problem ex- 
hibits t? dispersion. Numerical studies of chaotic incom- 
pressible flows typically show regions of chaotic and 
ordered behavior, therefore non-ergodicity. In order to 
study dispersion in chaotic regions, the usual procedure 
is to consider an initial distribution of points that is 
entirely contained in what seems to be a chaotic (and 
supposedly ergodic) region. From the above result, in 
the diffusive case it would not make any difference 
whether the points were initially placed in only one 
ergodic region, or spread out’ over several ergodic re- 
gions. The only importance lies in the fact that in the 
related non-diffusive problem, when the particles are 
placed initially in several ergodic regions, the nondif- 
fusive dispersion behaves like t* at large times. Then the 
Pe? regime is obtained in the diffusive problem. In this 
context we consider these results to be typical for a large 
class of laminar flows. 


Relative Time Scales of Lobe 
Transport Versus Transport by 
lecular Diffusion through Lobe 
Geometry 


So far, we have seen how, in the absence of molecular 
diffusion, lobe dynamics can be used to quantify the spreading 
of a passive tracer for finite times. Also, we have seen how 
ergodic theory can enable us to understand the asymptotic (in 
time) behavior of dispersion. In this section we show how the 
lobes can be used to determine a time scale over which lobe 
transport dominates molecular diffusion. Our knowledge of 
the dynamics of fluid particles associated with lobes suggests 
the following criterion that might explain the time scale over 
which lobe transport dominates molecular diffusion. 


The time scale, T,, for a tracer to diffuse across a distance of 
the order of a turnstile width should be long compared to the 
time it would take for a lobe to be mapped across the boundary 
of a region, i. e. T. 


Hence, we have 


(a(e))? 
Ty = D 


where d(e) is the maximum width of a tumstile lobe. Using the 
o(€) approximation to the tumstile width given by the Mel- 
nikov function, we have 





@ och(-- WA, 
Ty . (e 7 sech(> 7 )cosh( )*) 


ay @ 
T TD 





According to our criterion, lobe transport will dominate 
diffusive transport provided T, >> T. We check this for certain 
parameter values. With D=5.0x = and A= 0.1 we have the 
following three cases: 


1. @=0.62=0.1 implies Tp ~ 2007, 
2. @=0.24,e=0.1 implies Tp = 3007, 
3. @=0.6,e=0.01 implies Tp = 27, 

Thus, in the first two cases we would expect lobe transport 
to dominate molecular diffusion for about 200 and 300 peri- 
ods, respectively, and in the last case only for about 2 periods. 
This can be check by considering the generalized Langevin 
equations, which describe the influence of molecular diffusiv- 
ity: 


x= wv +n), (9) 


i= +C(0, (10) 


where 7(t),C(«) (diffusion terms) are random variables with a 
Gaussian probability distribution , such that their correlations 
are: 


<n(on(?) > = <C(NCG()>=2D&t- 7), (11) 
<n(OC(r’)> = 0. (12) 


In Figure 8 we compare the results of the lobe dynamics 
calculations (solid lines), which were obtained by numerically 
implementing the procedure described in Section 5, and the 
calculations which consider molecular diffusion(dashed 
lines), which were obtained by integrating (10). From this 
figure we can see that our criterion works reasonably well. 
Thus we see that geometrical and dynamical features of the 
lobes have lead us to a physical argument that gives us a time 
scale on which molecular diffusion can be neglected in the 
process of inter-cell transport. 


“Patchiness” in Fluid Flows 
and Distribution Functions 


Finally, we want to describe a new phenomenon that may 
prove to be a useful way of describing flows in a variety of 
geophysical settings. Pasmanter [1988], [1991] has studied 
mechanisms that give rise to the variability of dispersion 
processes in the ocean. A particularly important phenomenon 
to which he referred is known as “patchiness”, i. e. , a situation 


where parts of a distribution of passive tracer may disperse at 
different speeds compared to its surroundings. We want to 
show that the mathematical framework developed in Mezic’ 
and Wiggins [1994c] can be useful for studying this phenome- 
non. We will illustrate this with the cellular flow example. 

In the cell-to-cell transport problem the time average 
along particle trajectories of the x-component of velocity is a 
useful quantity. In order to illustrate the “patchiness” effect we 
numerically determine the distribution of average x-compo- 
nents of velocity for fluid particle trajectories in one cell. In 
Figure 9 we plot contours corresponding to initial points of 
trajectories that have the same average x-velocity which we 
obtained by numerically computing these time averages for a 
uniform grid of 1600 points for a time interval of 5000 periods. 
Regions of nonzero average x velocity correspond to “accel- 
erator modes”, i. e., regions of initial conditions corresponding 
to trajectories with nonzero average velocity in the x direction, 
and these are the points that participate in cell-to-cell transport. 
In Figure 10 we plot a histogram of this data, i. e. , we plot the 
number of points corresponding to a given “bin width” of 
average x velocity, where the bin width is 0. 002. “Patches” 
are mathematically defined as sets of positive measure in the 
flow on which the x component of velocity has a constant 
time-average along fluid particle trajectories. In Mezic’ and 
Wiggins [1994c] a distribution function is derived that allows 
one to deduce certain properties related to the patchiness. In 
particular, the “patches” correspond to points of discontinuity 
of a certain distribution function. An additional result obtained 
through study of this distribution function tells us that there 
are either a finite number of patches, or there is a countable 
number of them, but also countably many of them are of 
smaller measure than any prescribed number. In Figure 9, the 
“patches” are represented by the darkest and brightest spots. 
We could also define “patchiness” as the regime in which more 
than one “patch” exists. It then turns out that non-ergodicity 
of the flow is a necessary condition for “patchiness”. More 
precisely, it can be shown that if the flow is ergodic, there is 
only one “patch”, represented as one and only point of discon- 
tinuity of the distribution function derived in Mezic’ and 
Wiggins [1994c]. Figure 10 represents the “generalized de- 
rivative” of this distribution function, that is, the peaks corre- 
spond to different “patches”. From these results one can 
deduce a criteria for determining when a flow will be “patchy” 
If the flow is such that it has an ergodic component of positive 
measure, and if the initial distribution of passive tracer is such 
that the measure of that component is also of positive measure, 
then the flow will evolve into a “patchy” flow. 


Summary 


In this article we have taken a specific flow and shown 
how an arsenal of newly developed techniques from dynami- 
cal systems theory can be brought to bear on the study of the 


One 1995 13 





Lagrangian transport properties of this flow. This dynamical 
systems approach shows much promise in dealing with the 
types of Lagrangian transport issues that arise in geophysical 
fluid dynamics. Much of the promise lies in the fact that the 
mathematical framework of dynamical systems theory is re- 
markably similar to the experimental framework of modern 
flow visualization (from, e. g. , wacer data or satellite obser- 
vations) in that each is concerned with the motion in time over 
large regions of space and the role that lower dimensional 
organized “structures”, or geometrical features, in space play 
in governing the motion. Dynamical systems theory provides 
a mathematical framework for quantifying and describing the 
effects of these low dimensional “organized structures” on 
Lagrangian transport. Consequently, this type of mathematical 
approach should prove ideal for analyzing transport and mix- 
ing processes associated with the large scale, organized mo- 
tions observed in geophysical flows. 

The mathematical development of dynamical systems 
theory is also being influenced by studies of Lagrangian 
transport problems in geophysical fluid dynamics. The flow 
Studied in detail in this article was two-dimensional and time- 
periodic. Clearly dynamical systems techniques for three-di- 
mensional time-dependent flows must be developed if these 
methods are to have large impact on Lagrangian transport 
issues. In this setting, a variety of new mathematical issues 
must be addressed. The beginnings of the development of this 
theory can be found iin Beigie et al. [1991], [1992], [1994], 
Mezic’ and Wiggins [1994], and Wiggins [1992]. 

Finally, we mention a possibly revolutionary new devel- 
opment for Lagrangian transport problems. When the dynami- 
cal systems techniques are developed so that they can be 
coupled with modern computational methods for computing 
velocity fields the applicability of the methods will be greatly 
extended, as well as be more accessible to geophysical fluid 
dynamicists. Over the past 15 years computational fluid dy- 
namics has developed into a subject in its own right. Now we 
have accurate algorithms for solving the Navier Stokes equa- 
tions in a variety of physically important settings. However, 
many problems related to mixing and transport begin once this 
step has been accomplished. That is, first a solution to the 
Navier Stokes equations must be obtained in order to study the 
transport and mixing properties associated with that velocity 
field. In the vast majority of situations arising in applications, 
this solution can only be obtained numerically. Thus, we only 
have a numerical representation of the velocity field. Never- 
theless, many dynamical systems results do not care about the 
specific analytical form of the dynamical system under con- 
sideration. Rather, they require that only certain geometrical 
features be present. For example, the existence of stable and 
unstable manifolds of some invariant set requires only the 
existence of a hyperbolic invariant set and the existence of 
Smale horseshoe type chaos requires only the transverse inter- 
section of the stable and unstable manifolds of a hyperbolic 
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periodic orbit. Currently our research is focused on developing 
the theory and algorithms needed to perform dynamical sys- 
tems type analysis with numerical representations of vector 
fields. 
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However, the line design, line breakup, and line strength are discernable. On the cover, figures 8, 9, and 10 of this article 


are shown with color distinctions. 


1. Introduction 


The difficulties involved in collecting data on oceanic 
motion are evident and the inherent uncontrollability of the 
ocean for experimentation presents a formidable scientific 
challenge to those seeking to understand it. The data that is 
collected from the tracking of strategically placed floats fits 
very well with the Lagrangian perspective on fluid mechan- 
ics, wherein motion is viewed through the following of fluid 
particles, as it is assumed that the paths of these floats 
reflect those of fluid particles. These Lagrangian trajecto- 
ries are, from the dynamical systems viewpoint, the orbits 
generated by a vector field that is precisely the Eulerian 
velocity field of the fluid flow. This perspective thus offers 
a natural mathematical context through which to study 
questions about oceanographic flows. 

However, this form of oceanographic data collection 
through the tracking of floats suggests a deeper reason for 
the suitability of the dynamical systems perspective. Since 
only a small number of SOFAR (or RAFOS) floats can be 
placed in areas of the ocean under study, and they can only 
be tracked for a certain period of time, usually a year or less, 
the resulting particle paths must be viewed as a sampling of 
the corresponding Lagrangian trajectories. The pictures of 


these float paths are put together to form what are commonly 
known as “spaghetti” plots, see Figure 1 foranexample taken 
from reference 4. It is clear that we can only rely on general 
qualitative features of such plots for an assessmentof the fluid 
motion. One cannot hope to give an explanation of detailed 
properties of these paths. Indeed, such an explanation would 
not be meaningful as any individual path is the inevitable 
victim of small-scale, temporary effects. However, the large- 
scale features of the paths that are present in a number of 
samplings should beexplained. 

The modern theory of dynamical systems, with its 
emphasis on qualitative properties and characteristic geo- 
metric features of trajectories, provides a fitting perspective 
on the issues raised by the study of oceanographic flows. 
The underlying tenet of geometric dynamical systems the- 
ory is that certain basic geometric structures can be isolated 
which organize the overall structure of the flow. In particu- 
lar, certain specific trajectories, if they can be located, will 
demarcate different regions of qualitatively distinctive tra- 
jectory motion. Specifically, stagnation points of the flow 
and their attendant stable and unstable manifolds (sets of 
trajectories approaching the stagnation point in, respec- 
tively, forward and backward time) play definitive rdles in 
organizing the flow. 


One 1995 17 





From the dynamicist’s perspective, the oceanographic 
flow fields are presented to us on three levels of increasing 
accuracy, and complexity. The simplest level of model is of a 
velocity field that is postulated as an accurate model of the 
phenomenon under investigation. This will be formed, com- 
monly, in terms of an unstable, or neutral, mode superimposed 
on a base flow. Such a flow leads to an analytic expression for 
an ordinary differential equation whose analysis is, in princi- 
ple, possible by existing techniques and only restricted by its 
complexity, in particular its dimension. The next level of 
model is the numerical model which consists of a velocity field 
given as a numerical database that is constructed by solving 
the relevant partial differential equation on an appropriate 
domain. A dynamical systems approach to such a velocity field 
is severely impeded by the fact that all the developed tech- 
niques are based on the idea that the velocity field is supplied 
to the mathematician analytically. A new perspective is thus 
called for and the development of techniques for handling 
numerically generated velocity fields is a significant task that 
remains to be adequately addressed. In particular, the numeri- 
cal determination of the relevant structures, such as stable and 
unstable manifolds, for numerical velocity field data sets is 
required. A strong case has been made by Ridderinkhof and 
Zimmerman”, for the relevance of these dynamical structures 
in a numerical model of the Wadden Sea. The final level is the 
experimental data set available for the situation under inves- 
tigation. These data sets, as discussed above, are usually based 
on the tracking of floats and, while insufficient to be viewed 
as a complete model, must always be used as the final arbiter 
of the success of a given piece of analysis. 

In this review, we shall consider models at the first level. 
The formulation of such models as base flow plus perturbing 
mode can offer great insight into charcateristic flows that may 
truly exist as the nonlinear equilibration of such “neutral” or 





Figure 1. 


Spaghetti diagram of RAFOS floats in the northern Atlantic Ocean 
off the eastern United States from 1984-1985. Figure taken from 
Bower and Rossby* 
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“unstable” directions from the base flow. Explicit anaiytic 
expressions, or even approximations, for such fully nonlinear 
Eulerian flow fields are comparatively rare, and these “linear” 
approximations offer, with the existing technology, the best 
way to make meaningful, analytic predictions about the nature 
of Lagrangian trajectory motion in the full flow field. How- 
ever, a goal of the mathematical development motivated by 
oceanographic problems must be the construction of better 
analytic models that incorporate real nonlinear effects while 
remaining faithful to the underlying physics. 

Our focus will be on heteroclinic orbits, which are trajec- 
tories joining different stagnation points, and are often known 
as separatrices on account of their function in separating 
different flow behaviors. We favor here the more modern term 
of “heteroclinic orbit” because the main point of the message 
is that, under perturbation, these structures no longer perform 
this function of separation. Nevertheless, there remain, typi- 
cally, heteroclinic orbits after perturbation. We shall present 
three examples that illustrate how an elucidation of the struc- 
ture of these orbits, in particular for small perturbations of the 
base flow, can reveal the mechanisms in the flow for the 
transport of particles. 

The first example concerns the settling of a particle in a 
cellular flow field. The particle of interest here is not, in fact, 
a fluid particle but an aerosol particle, and a key point in the 
analysis is the small inertial effect experienced by such a 
particle. the question of interest is whether the particle under 
consideration can be trapped in the cell or continue to settle 
through the flow indefinitely. 

An interesting feature of jets in the ocean, such as the Gulf 
Stream, is the spawning of rings that form a vortex motion, 
and this is the motivation for our second example. 

These rings carry, for instance, cold water into the region 
of warm water to the south of the Gulf Stream. They tend to 
travel in the reverse direction from the jet and survive for a 
significant amount of time, often unul hitting the continental 
shelf. The formation of these rings is not well understood, see 
Pratt et al.”° and many of their features remain to be explained. 
Kirwan et al.'*!® formulated a model for such rings in terms 
of rotating modons. These are special solutions of the qua- 
sigeostrophic equations and their mixing properties will be 
discussed in section 3. 

The third example concerns the transport of water within, 
and across, an oceanic jet. The general mechanisms for mixing 
across a jet, such as the Gulf Stream, are of fundamental 
importance. Various kinematic models have been formulated 
for the clarification of this transport, but, recenuy, models of 
velocity fields that are dynamically consistent have been de- 
rived by Pratt et al.2! We discuss in section 4 some of the issues 
involved in the transport of fluid across such a jet. 

In each of these examples, there is an “unperturbed” flow 
field that is a one-degree-of-freedom Hamiltonian system. The 
distinguished orbits that connect stagnation points to each 








Figure 2. 


A diagrammatic representation of the streamlines of Stommel’s 
model. The protected pods can be seen in each of the four cells 
depicted. Trajectories in these pods correspond to particles that 
are forever suspended, while the particles outside will settle. 
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other separate regions of contrasting Lagrangian motion from 
each other. The perturbations supplied by effects unaccounted 
for in the base flow will tend to break these separating orbits 
and offer the possibility of transport between regions that 
were, in the unperturbed case, characterized by different mo- 
tions. We shall see, in each of the examples discussed below, 
the implications of the perturbations on these distinguished 
orbits. 


2. Settling in a Cellular Flow Field 


The problem of whether particles of non-fluid matter 
settle in the ocean under the influence of gravity is obviously 
of great importance. A cellular flow field is a model of Lang- 
muir circulation and we consider here the issue as to whether 
aerosol particles will tend to settle through such a flow field 
or be trapped and recirculate within that cell. A key point is the 
réle of the inertia of the particle in question. The inertial effects 
can provide, as will be seen below, an effective dissipation of 
the Hamiltonian that arises as the streamfunction of a 2-dimen- 
sional incompressible flow. 

The settling of an aerosol particle through a cellular flow 
field was studied by Stommel”’. Stommel’s analysis neglected 
the inertia of the particle, and thus the particle was effectively 
treated as a fluid particle settling through the flow field under 
the combined influences of gravity and the surrounding fluid. 


He showed that, in this scenario, both trapping and settling 
were possible. Indeed, the ratio of trapped particles to settling 
particles varied with the characteristic terminal settling veloc- 
ity due to gravity. In physical space, when the influence of 
inertia is omitted, heteroclinic orbits to the stagnation points 
on vertical invariant lines between cells separate the regions 
from which particles will settle from those in which particles 
will be trapped, as shown in Figure 2. 

Indeed, these heteroclinics form invariant, protected 
pods. Under the periodic velocity field, the Lagrangian trajec- 
tories outside a pod will settle through all of the cells they 
enter, whereas the Lagrangian trajectories inside a pod will be 
trapped forever. 

Despite the apparent possibility of many particles being 
suspended in the cells, Maxey and Corrsin'*-'* observed, based 
on numerical experiments, thai inertial particles had a greater 
tendency than fluid particles to settle. The true problem with 
inertia included is, however, a perturbation of that considered 
by Stommel, and we must consider whether any analogue of 
the protected pod, which might guaranteee the suspension of 
particles, exists in the inertial problem. 

If the inertia of the particles moving through the cellular 
flow field is small, but nonzero, the velocity field for these 
particles becomes a singularly perturbed version of that for 
fluid particles acted upon by the cellular flow field and gravity. 
The analysis and interpretation of such singularly perturbed 
problems contains many pitfalls. Fortunately, a technology 
exists which can render a complete geometric description of 
the trajectories of such perturbed systems, and it can be in- 
voked here to give a clear and unambiguous answer as to 
whether, under the influence of inertia, particles will prefer to 
settle or not. 

This set of techniques is based on a striking collection of 
theorems due to Fenichel.’*:'* These results give conditions 
under which, for the perturbed system, there exists an invariant 
manifold (a smooth surface of solutions) that carries a vector 
field which is a small (regular) perturbation of the formal limit 
equation. 

If we let y(t) denote the position of the center of the 
particle in space (y=(y,,y>)) and take the downward direction 
to be the direction of positive y, (the opposite of the conven- 
tion in reference 11) then the equations of the following model 
govern the motion of an 2erosol particle in a cellular flow field 


y=) 


y=", 


€v, =— v, +Sin y; COS y> 
€¥7+—V2—W—COSy, SiNy, (I) 


Here, v(t) is the velocity of the particle and W is the 
settling velocity scaled so that 0 < W< 1. The small parameter 
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Here, v(t) is the velocity of the particle and W is the 
settling velocity scaled so that 0 < W< 1. The small parameter 
€ > 0 is the Stokes number and measures the inertial response 
time of the particle to the medium. The case ¢ = 0 renders the 
formal limit in which the particles have no inertia and behave 
merely as fluid particles. Fenichel’s work implies that there 
exists, for system (1) with small inertia, an invariant manifold 
of the type discussed above. This invariant manifold, known 
as the slow manifold, is given by a small perturbation of the 
algebraic conditions found by formally setting € = 0 in the last 
two equations of (1). An application of Fenichel’s further 
results combined with some simple dynamical systems argu- 
ments show that, in this case, the slow manifold is actually 
globally attracting. Thus the asymptotic motion of all La- 
grangian trajectories will be determined by the vector field on 
the slow manifold. 

A key point to notice is that the slow manifold is para- 
metrized by the physical space variables (y,,y,). To analyze 
the true behavior on the manifold, it is necessary to calculate 
the O(€) terms in the perturbed vector field. This is achieved 
by using the invariance property of the manifold and the full 
equation itself. The resulting system on the slow manifold is 


y; =Sin y, cos y, —€ sin y, (W sin y, + cos y,) + O (e”) 


(2) 


¥2 =— W—cos y, sin y, +€ cos y, (W cos y, +sin y,) + O(e*). 


As commented above, the long-time behavior of aerosol 
particle trajectories is determined completely by the trajecto- 
ries of this vector field. Moreover, the interpretation of these’ 
trajectories is considerably facilitated by the fact that the 
dependent variables in (2) are exactly the space variables. 
Setting ¢ = 0 in (2), we recover the velocity field for fluid 
particles, the trajectories of which are sketched in Figure 2. 
We are now ready to address the key question as to whether 
particles should be expected to settle under the influence of 
the vector field (2). The answer will obviously come from 
studying the fate of the pods present in the ¢ = 0 system under 
the perturbation supplied by the O(e) terms of (2). This is 
tantamount to asking how the heteroclinic orbit surrounding 
the pod in the interior of each cell will break under the 
influence of this perturbation (the vertical heteroclinic stays in 
place). A few simple sketches will convince the reader that if 
the unstable manifold falls outside the stable manifold then the 
process of settling is facilitated, whereas the reverse configu- 
ration would result in more trapping. The cc:ermination of the 
direction of breaking of heteroclinic orbits is the objective of 
Melnikov calculations,'°*° and the issue can be resolved by 
calculating the sign of a certain integral. This integral presents 
many technical difficulties, and we do not evaluate it directly, 
but itcan be shown analytically that its actual sign corresponds 
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to the unstable manifold falling outside the stable manifold. 
This indicates a strong enhancement of settling'’. 

In this way, we show that any small inertial effects in fact 
force particles to settle through a cellular flow field and make 
the probability of suspension almost non-existent. 

Figure 3 provides a sketch of the resulting flow of parti- 
cles with small nonzero inertia. 

It is also shown analytically in reference11that certain 
characteristic paths emerge in the slow manifold, along which 
the aerosol particles will, asymptotically, congregate and 
travel. This corroborates further numerical observations due 
to Maxey and Corrsin.'>! 

Although the underlying cellular flow field is a reasonable 
model in certain oceanographic circumstances, it is, of course, 
very idealized. Nevertheless, our work shows that the inclu- 
sion of inertial effects, if properly interpreted on a slow mani- 
fold, supplies an effective dissipation of the Hamiltonian 
streamfunction. Such effects are also seen in recent related 
work of Tio et al.”*, in which settling through a critical layer 
is studied. This may partly account for observations in which 
particulate matter in the ocean is seen to collect in certain 
well-defined areas”?**, despite the fact that it is under the 
influence of a time-dependent, incompressible, almost strati- 
fied flow, which one would otherwise expect to be strongly 
homogenizing. 


3. Rotating Modons 


In recent work, Mied, Kirwan and Lindemann'® have 
proposed rotating modons as models of rings that result from 





Figure 3. 


The trajectories of those particles lying on the slow manifold are 
depicted here. Other trajectories will asymptotically approach this 
configuration. Small inertia is included and has the effect of 
breaking the heteroclinic orbits in such a way as to force settling. 
Particle paths do not become trapped and those lying within the 
protected pods in the zero inertia case can leak out and settle. 




















Figure 4. 


Hamiltonian (effective streamfunction) of the feature model for 
82-B. Solid and dashed curves are positive and negative values, 
respectively. The contours have been normalized by aconvenient 
plotting scale. The dotted circles show model scale lengths for 
the upper and lower model layers. The inner circle radius is 120 
km. There are two types of critical points in this plot. The three 
diamonds depict stagnation points that are local extrema. There 
are two dots near the tips of the SP contour lines. These are saddle 
points. The SP contour is the value of teh Hamiltonian at the 
saddle points. 




















the pinching off of vortex-like regions from such jets as the 
Gulf Stream. These multipole vortex systems consisting of one 
or more cyclonic eddies moving around an anticyclonic ring 
often are seen in AVHRR imagery in the coastal/open ocean 
interface. The nonlinear feature model developed by Kirwan 
et al.'* accounts for many of the charcteristics seen in the 
imagery. These modons are specific solutions of the qua- 
sigeostrophic potential vorticity equations for a 2-layer model, 
in such a setting they were first introduced by Flierl et al.?. 
Modons are steady in a frame that is rotating at a fixed 
frequency. Bottom topography can be included, and this can 


have the effect of trapping the modon and preventing it from 
lateral movement, see’*. Of particular interest are dipolar 
vortices in which two regions of opposite rotation are matched 
together. It is suggested in reference 16 that such dipoles model 
pairs of warm and cold cores that are attached to each other, 
and can even become so attached as the result of a collision. 
Of interest in our study is the local mixing properties of such 
modons. For instance, is it possible that fluid from outside the 
modon can be mixed into it? Can it flow into the core? Or 
between the two complementary cores? If not, what external 
mechanisms might facilitate such mixing? 

The rotating modons are solutions of the quasigeostrophic 
equations of a prescribed form. They are constructed by find- 
ing solutions in different circular regions and matching at the 
boundaries. This is quite a complicated procedure that leads to 
formidable calculations for solutions of the resulting partial 
differential equations. 

In the paper by Kirwan et al. '*, each layer is divided into 
an inner (circular) modon region and an outer region. The radii 
of the regions in each layer will, in general, be different. The 
quasigeostrophic equations for the potential vorticity q; and 
the streamfunction w,; are assumed to hold in the regions of 
each layer, i = 1,2, with appropriate matching conditions at the 
boundary. Solutions of these equations which are steady in a 
rotating frame of frequency @ can be found under the ansatz 
that the effective streamfunction 





Figure 5. 


This vortex system has a rotation period of 30.6 inertial days. 
Positive values are depicted as solid curves and negative 
values as dashed curves. Two blobs are initialized on either 
side of a saddle point; the red blob is just inside the saddle 
point while the biue blob is just outside. The blob radii are 5km 
and there is a 2 km distance between the blobs. 























Figure 6. 


Flow associated to the single cusp jet model with a single neutral 
mode superimposed. When viewed in a reference frame moving 
with the speed of the neutral wave, the flow is time-independent 
and admits two rows of cnitical points at the steering lines. Fluid 
between the steering lines moves to the right as its speed is 
greater than the speed of the reference frame, and fluid 
outside the steering lines moves to the left. In each of the 
figures of this section, the heteroclinic orbits are colored red 
(blue) for the pieces that are generated as unstable (stable) 
manifolds. This is not so meaningful in the unperturbed case 
(Figures 6, 7 and 9) since, on the heteroclinic orbits the stable 
and unstable manifolds coinicide. the reamining trajectories 
are colored red. 



































or 

Yi=Vi- > 

and the potential vorticity are functionally related, with a 

different relation in each layer and region. The quasigeostro- 
phic equations are then 


J(¥,4) = Sj, 


where S; are sources at the boundaries of the regions to 
compensate for the jumps in potential vorticity. The relevant 
solutions are then found by Kirwan et al.'* from solving the 
partial differential equations coming from the functional rela- 
tions. 

In physical variables, namely r, 8 and t, the Lagrangian 
trajectories are those of a Hamiltonian system with time-de- 
pendent forcing. There are two such systems, one for each 
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layer, and the Lagrangian trajectories in each layer are inde- 
pendent of those in the other layer, although the corresponding 
streamfunctions are not. Within each layer the Lagrangian 
trajectories are the same for any depth. The trajectory equa- 
tions are then 


dr_ 1 OW; 


dt rae 


Ci) 
dt 


Vi 
or’ 


i 
r 


where the streamfunction y, _ w,(r,6 - wt). The streamfunction 
will be continuous across the boundaries of the various re- 
gions, but not necessarily smooth. Since the system (3) is a 
time dependent Hamiltonian system, its trajectories can be, 
ostensibly, chaotic. However, if viewed in a rotating frame, it 
is easily seen that the effective Hamiltonian Y; becomes a time 
independent Hamiltonian and after the appropriate changes of 
coordinates, (3) can be recast as a one-degree-of-freedom 
Hamiltonian system. We introduce the variables © = @ - wt and 
a scaled time t, so that (3) becomes 





Figure 7. 


Increasing the wavenumber of the neutral mode draws the 
steering lines closer to the centerline of the jet and the meander 
of the jet becomes more pronounced. The assymetry in the 
cat's eyes is due to the exponential form of the eigenfunctions. 
Here the neutral mode has wavenumber k = 8 and amplitude 
a=0.01. 
































Figure 8. 


Adding a second neutral mode to the jet model results in a time 
periodic vector field. Here the second mode is moving faster 
than the first mode with wavenumber k2 = 16 and amplitude € 
= 0.001. In the presence of this time periodic perturbation the 
heteroclinic orbits break up into transverse intersections of the 
stable and unstable manifolds for the associated Poincare 
map. The most significant mixing occurs on the inside of the 
jet, however, no mixing occurs across the body of the jet. The 
color coding (red-unstable, biue-stable manifolds) shows 
clearly the intermingling between stable and unstable mani- 
folds on the inner heteroclinic orbits. 




















do OY; 


dt or 


where the right hand side depends only on the phase space 
variables r and ©. 

This casting of the problem has rather striking implica- 
tions. If the modon is viewed in the rotating frame, then the 
trajectories must follow the contours of the effective Hamil- 
tonian and are, in particular, quite tame. We shall focus on the 
case of a dipolar modon. In recent, and ongoing, work Kirwan 
and the authors have discovered stagnation points in the flow 
field (3) for the upper layer. These correspond to periodic 


Lagrangian trajectories in physical space. These stagnation 
points are crucial in understanding the dynamics of (3). Indeed 
they are located near the edge of the modon boundary (radius 
of inner region) for the upper layer. The outer ring of hetero- 
clinic orbits (called a heteroclinic cycle) surrounds both the 
cyclonic and anti-cyclonic vortices (warm and cold cores) of 
the modon. The inner cycle surrounds only the cyclonic vor- 
tex. Figure 4 shows the qualitative nature of the phase portrait 
of (3). The figure is based on data for a model of 82-B, using 
hydrographic and AVHRR data. This data is well documented 
in Kenelly et al.'’ The two stagnation points are located in the 
north-east and north-west sides of the modon structure. The 
heteroclinic orbits are contours of the effective Hamiltonian 
joining the stagnation points. These are not cpatured exactly 
in the Figure 4, but reasonable approximations can be seen in 
the contours marked SP, it is left to the readers imagination to 
see how these are connected to the saddle stagnation points to 
form the heteroclinic orbits. Of particular significance in this 
analysis is that the heteroclinic cycles provide a barrier to the 
mixing of fluid in certain regions. The outer homoclinic cycle 
prevents any outside fluid from infiltrating the modon struc- 
ture. This is not obvious from pictures in the physical space as 
this homoclinic orbit will not be circular and hence will change 
position as the modon rotates. Nevertheless, if a fluid particle 
is outside the rotated version of this cycle at an appropriate 
time, the modon flow will never carry it through the modon. 
In some sense, this cycle defines the modon as it supplies the 
modon boundary. The inner heteroclinic cycle provides the 
barrier between the two main vortices of the modon. It should 
be noted that the contouring of the effective streamfunction 
(Hamiltonian) has uncovered a third vortex patch which was 
somewhat unexpected. It is, however, viewed as less signifi- 
cant than the two primary vortices. 

A numerical experiment, performed by B.Lipphart (Old 
Dominion University), gives clear corroboration of the protec- 
tion of the modon by the outer heteroclinic cycle. In this 
numerical experiment, shown in Figure 5, two blobs of initial 
data are placed in locations near the northeasterly stagnation 
point. The outer blob (blue) is placed just outside the outer 
heteroclinic cycle at the stagnation point, and the inner blob 
(red) just inside. The resulting particle paths are then strobed 
at six, equally spaced, time intervals in physical space ending 
after one period. The contour lines sketched in Figure 5 are 
those of the original streamfunction (y, and not Y;) and their 
presence is used to facilitate the viewing of the modon. The 
different fates of the two blobs can be seen quite clearly from 
the numerical plots. Indeed, the blue particle paths are kept out 
of the modon structure, while the red ones are sucked into the 
region between the vortices. This difference is especially 
striking when one considers that the initial separation of the 
blobs is not that great, only about 2km. in physical space. 

In the real oceanographic systems, one cannot expect the 
modons to be totally protected from the outside fluid. Indeed, 
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the outside fluid will be sucked through the center of the 
modon and, to some extent, into the edges of the vortices. The 
question then arises as to what perturbing mechanisms might 
effect such mixing of fluid through the modon structure. This 
is the subject of our ongoing work in collaboration with 
Kirwan, but some observations can be made as to why and 
how we anticipate this being realized. The outside ocean is a 
ready supplier of perturbations that will manifest themselves 
as time varying influences on the modon structure. It is a 
well-known fact in dynamical systems that for area-preserving 
systems, homoclinic orbits for the associated Poincaré map 
will persist under perturbation. This can be applied here but, 
since the perturbed system will be time-dependent, has the 
consequence that the stable and unstable manifolds of a fixed 
point, close to the original stagnation point, will intersect. 
Since these are now invariant manifolds for a map, their 
intersection produces a quite different scenario. Indeed, the 
manifolds for the perturbed system should intersect trans- 
versely, which we expect as a consequence of a symmetry 
argument similar to that noted in the following section, and 
chaotic mixing of fluid particles between the regions that were 
previously protected from each other is now likely,”'. It may 
be commented here that the conservation of potential vorticity 
will prevent chaotic mixing, as observed recently by Brown 
and Samelson*. However, in the unperturbed case, the poten- 





Figure 9. 


Modeling the dynamics at a steering line where 

U'(y) = 0. Here the steering line occurs at a point where the 
base velocity profile attains a minimum (U” (0) 0). Because of 
the scaling in this analysis the dynamics actually occur over a 
smaller range of y than in the previous figures. 




















24 Naval Research Reviews 


tial vorticity and effective streamfunction are linearly depend- 
ent, that is indeed how they were found for the modon, and 
thus the analysis of reference 5 does not apply in this case. 
Even though we expect mixing through the modon to occur 
under such time-dependent perturbations, the KAM Theorem, 
see reference | for instance, can be applied to show that, under 
small perturbations, some of the orbits surrounding the vortex 
regions will survive as invariant tori. This will protect the vortices 
themselves from being sucked into the mixing process. Thus, we 
expect for the flow that fluid will be mixed through the modon, 
that is to say between the two vortex cores, but not into the cores 
themselves. 

A key point here is how to model the time dependent 
perturbation, whether by a superimposed mode or a time 
variation of parameters in the problem, such as the radii of the 
regions in one of the layers. Issues of dynamical consistency 
and physical relevance are of central importance. However, 
we expect that mixing of fluid will occur through the modon 
for almost all such perturbations. Given the formulation of an 
appropriate model, mixing rates can even be calculated from 
the Melnikov function, as formulated by Rom-Kedar and 
Wiggins”*>', and for some more recent theory”, and the 
different mechanisms available can be compared for the effi- 
cacy in the process of mixing through the modon. 


4. Mixing across Jets 


The mechanism for the transport of fluid particles across 
such oceanic jets as the Gulf Stream is not well understood. 
Near the surface the presence of high potential vorticity gra- 
dients at the center of the jet appear to act as barriers to this 
transport. Indeed, this has been experimentally corroborated 
by Behringer, et al.”. Nevertheless, it appears to occur either 
by ring detachment or straightforward advection. 

Many models have been developed with respect to which 
these questions can be posed. Bower’ and Samelson”® studied 
kinematic models based on the meander structure of the Gulf 
Stream evident from spaghetti plots, such as shown in Figure 
1. More recently, dynamically consistent models were formu- 
lated by del-Castillo-Negrete and Morrison® and Pratt, Lozier 
and Beliakova”'. The first group used the Bickley jet (sech’y 
profile) as the base flow coupled with superimposed neutral 
modes. Pratt et al used the cusped jets derived from a 12 layer 
model in the paper by Pratt and Stern”. These cusped jets have 
piecewise constant potential vorticity, and the center of the jet 
is characterized by a discontinuity in the potential vorticity. 
Since these cusped jets are linearly stable, it is reasonable to 
model perturbations in the full flow field by superimposing on 
the base flow selected neutral modes. It is assumed that two 
waves are used and the first is a larger perturbation than the 
second. The first neutral mode introduces a characteristic 
speed and the steering lines occur where this matches the speed 
of the fluid particles in the base jet. As pointed out by both 








Figure 10. 


Adding a second neutral mode with a speed slightly greater than 
the first mode results in significant mixing across the region of the 
steering line. This is evident from the strong breakup of the stable 
(blue) and unstable (red) manifolds of the associated Poincare 
map. Again the intermingling is shown by the color coding, see 
Figure 7. 




















groups, in the regions near the steering lines chaotic trajecto- 
ries Can appear and, more importantly, advective mixing be 
facilitated. Such chaotic structures will be expected to appear 
when a second wave is superimposed on the flow field. If its 
wavenumber and frequency are unmatched to those of the 
primary wave, the resulting forcing on the velocity field is time 
dependent and such a scenario is known to produce chaotic 
motion. 

A comment about potential vorticity is in order here, see 
the beok by Pedlosky’® for the relevant definitions and discus- 
sion. As observed recently in the paper by Brown and Samel- 
son’, a two-dimensional incompressible flow with 
time-periodic forcing, which preserves potential vorticity, 
cannot be chaotic if the gradients of the potential vorticity and 
the streamfunction are independent. Indeed, the Lagrangian 
trajectories are constrained in this case to 2-dimensional sur- 
faces, determined by the potential vorticity contours, and will 
thus not be chaotic. As commented by del-Castillo-Negrete 
and Morrison®, this apparent contradiction is resolved by the 
fact that flows under consideration are not exactly potential 
vorticity conserving. However, from the mathematical view- 
point, this breaking of the potential vorticity conservation 
arises from the approximations made in implementing the 


model. These approximations are made in assuming that the 
flow field is given as a base flow with two neutral modes 
added. The full nonlinear field, which is only modelled by this 
linear superposition, is potential vorticity conserving, but this 
may not hold for the flow field that is actually used. Moreover, 
potential vorticity can be broken by the numerical approxima- 
tions made in calculating the Lagrangian trajectories. Of 
course, in the real ocean potential vorticity is not conserved 
and we do expect orbits heteroclinic to stagnation points to 
break up in such a way as to produce chaotic motion. It would 
be preferable, however, if this effect was modelled directly 
rather than as a by-product of the approximations made. In 
other words, we would hope to be able to produce these effects 
of chaotic mixing by producing models that more accurately 
reflect the actual ways in which potential vorticity is broken, 
rather than through less accurate ways of modelling systems 
in which it is conserved. Nevertheless, the presence of chaotic 
mixing is credible whenever such heteroclinic orbits appear 
and we believe that it is a significant mechanism. 

The models under consideration are effectively two-di- 
mensional, and the effects of vertical motion may be quite 
significant. It has been suggested that mixing across a jet may 
occur more readily at deeper levels, see Lozier and Riser’. 
Combining this feature with the possibility of vertical mouion 
offers some very tantalizing possibilities. The effects of depth 
can be studied within the stratified models discussed above as 
the depth does occur as a parameter within the flow field. 
Questions can then be asked about the dependence of advec- 
tive mixing on depth, for instance: is mixing facilitated at 
lower depths? The picture one should keep in mind is that 
described by Pratt, Lozier and Beliakova”' in which the base 
jet is bracketed by two steering lines determined by the pri- 
mary wave. On the surface, these steering lines are well 
separated, but as one looks deeper, the steering lines move 
toward each other unul they converge and cancel each other, 
at a certain depth. Below that depth the fluid particles of the 
base jet all move slower than the primary wave. 

To study the dynamics of the Lagrangian trajectories in 
the vicinity of steering lines we began with the analytical 
model for a cusped jet velocity profile formulated by Pratt, 
Lozier and Beliakova”'. The flow is taken to be two dimen- 
sional and is governed by a streamfunction y(x,y,t). The 
Lagrangian particle motions satisfy the equation 


de __ oy 


dt ay 


dy _dy 
dt ox 


The base flow of the jet is a shear flow with velocity profile 


Uy) =e™. 
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W(x.) = a@P(y) cos k(x — ct), 
where the eigenfunction (y) is given by 


My) =e*™. 


The wavenumber k, the wave speed c, and the cross-stream 
decay rate « are related by 


| 


K=V1+EF, c+1-« 


In a frame moving with speed c, the Lagrangian trajecto- 
ries are solutions of the equations 


4 =U y)-c—a 0) cos Kf, 


a , 
te ak® (y) sin kC, 


where we have introduced the variable C = x - ct. As shown in 
21 this system will have two steering lines, for appropriate 
values of c, and the structure of the trajectories near these lines 
will be of a row of cat’s eyes, familiar from the Stuart vortex 
flow, see Figure 6. Increasing the speed c of the superimposed 
wave draws the steering lines closer to the centerline of the jet. 
The exponential form of ®(y) introduces asymmetry into the 
shape of the periodic orbits and the meandering shape of the 
jet is more pronounced, as seen in Figure 7. 

The break-up of these heteroclinic trajectories is forced 
by a second neutral wave being superimposed on the stream- 
function. This secondary wave is of the same form as the first 
wave, 


W2 (x,y,t) = €®4(y) cos ky (x - cf), 


but the wave speed c, is chosen different from the speed c of 
the first mode. The wavenumber k, and speed c, are related as 
in (5), and the amplitude ¢ is taken to be small relative to a. 
The Lagrangian trajectories are now governed by a time-de- 
pendent vector field where the order € time dependency is 
periodic with angular frequency @ = k,(c,-c). The vector field 
is of the form 


&-uy) —c—a®’ (y) cos kt — €®’, (y) cos (KC — wr), 


(7) 


D =—ak® (y) sin kG —€ k > (y) sin (kyl - wr). 


Since (7) is a time-dependent Hamiltonian system for € 
0, the secondary wave is expected to break the heteroclinic 
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orbits on either side of each set of cat’s eyes. This can be 
verified by calculating the Melnikov function associated with 
the problem. This function gives the leading order term in the 
signed distance between the stable and unstable manifolds of 
two unstable periodic solutions which are created by the 
perturbation in the vicinities of the two unperturbed stagnation 
points. For the vector field (7) the Melnikov function can be 
written in the form 


M (t,) ” f- {Wo - Var¥2}, . ° (8) 


where x"(t-tg) is the unperturbed heteroclinic orbit connecting 
the stagnation points for ¢ = 0, and {.,.} denotes the Poisson 
bracket in the C and y variables. The explicit form of x*(t-tg) 
seems hopeless to determine, but it is in fact not needed for 
finding zeros of the Melnikov function. It turns out that, using 
the symmetry of x"(t-t,) with respect to the fixed € line through 
the center of the unperturbed problem, one can argue that M(t) 
admits isolated zeros. A similar reasoning yields that these 
zeros are transverse, therefore the heteroclinic connection 
between the unperturbed stagnation points breaks up, but 
individual transverse heteroclinic solutions between nearby 
unstable periodic motions survive. 

This break-up can be observed numerically through the 
distinctive looping caused by transverse intersections and the 
associated heteroclinic tangle. Figure 8 shows the result of 
perturbing the flow shown in Figure 7 for the case when the 
wavenumber k, of the secondary mode is twice the wavenum- 
ber of the first neutral mode. An interesting feature here is that 
it is significantly harder to detect the presence of transverse 
intersections along the outside edge of the jet. This is similar 
to a feature noted by Samelson”® for the kinematic meander 
model. The interesctions on the outer edge seem to occur at a 
much higher order, and it appears that the outside heteroclinics 
do not break at all. Whereas, when the steering lines draw 
closer to one another, the heteroclinics in the interior of the jet 
are easily detected. In agreement with the analysis for the 
derivative of the Melnikov function, we have observed nu- 
merically that the strong heteroclinic breakup can be made to 
occur on the opposite side of the loop by reversing the sign of 
the angular frequency @. This implies that the regions of 
significant mixing may depend on the relative speeds of the 
two superimposed neutral waves, since the sign of w depends 
on the difference of the wave speeds c and c,. In particular, if 
the secondary mode is faster than the primary neutral mode, 
then the resulting mixing is more intense within the body of 
the jet, whereas the boundary of the body of the jet, viewed as 
being the two steering lines, is virtually intact. As a result, 
fluids particles are kept outside of the jet by the outer stable 
and unstable manifolds which form barriers to mixing. At the 
same time, if the wave speed of the secondary mode is smaller 
than that of the primary mode, this boundary breaks up and 





being the two steering lines, is virtually intact. As a result, 
fluids particles are kept outside of the jet by the outer stable 
and unstable manifolds which form barriers to mixing. At the 
same time, if the wave speed of the secondary mode is smaller 
than that of the primary mode, this boundary breaks up and 
Lagrangian particle motions can penetrate into the jet via 
chaotic mixing. However, even in this second case the inside 
of the jet is bounded by the virtually unbroken inner separa- 
trices and serves as a barrier to transport through the core of 
the jet. 

As steering lines approach each other, which occurs as 
one goes deeper below the surface, it appears that mixing, both 
within and across the jet, should be facilitated. Our group is 
studying this in the context of the cusped jets. For the jet with 
a single cusp, the steering lines approach each other, but the 
limit is quite degenerate as they will converge exactly at the 
cusp. The jet clearly becomes very thin, but indications are that 
it remains a barrier, as seen in Figure 8. 

In the Bickley jet case, ° the phenomenon of separatrix 
reconnection occurs. In this scenario, heteroclinic orbits exist 
that stretch across the jet between the two rows of cat’s eyes 
that lie either side of the jet centerline. This would appear to 
be a level at which mixing is strongly facilitated and can be 
realized across the entire jet. However, as pointed out by Pratt 
(private communication) this may be hard to realize in an 
actual situation because of the strong potential vorticity gradi- 
ents. Indeed, this is corroborated by the single cusped jet case. 

The phenomenon of convergent steering lines can also be 
studied in the context of the double cusped jet of Pratt et al.”. 
Here they will converge as the speed of the primary wave is 
reduced to match the speed at the center of the jet, which lies 
in a third region of constant potential vorticity between the two 
outer regions. In the case of a smooth jet profile, the derivative 
U’(y) is necessarily zero at the speed at which the two steering 
lines converge. We have performed a local analysis to study 
the possible dynamics in such a situation. For simplicity, we 
assume that the converging steering lines occur at y=0 and that 
both {y) and ®’(y) are non-zero at y=0. We introduce the 
following change of variables: 


C=x-ct, n=a-“y, t=a%. 
Retaining just the leading order terms of the vector field results 


in the following scaled equations, 


a =} U” (0) y- — a” & (0) cos kt, 


(9) 


AD = — k@ (0) sin K—€ a ky ©, (0) sin (ky C—a on). 


Ate=0 there are homoclinic loops separated by a distance 
kC = 2x, and heteroclinic orbits form separatrices on only one 


side of the jet. Figure 9 shows the resulting flow with one 
neutral mode superimposed and the parameters set to U”(0)=2, 
(0)=1, ©(0)=1, and a= 0.01. A significant issue is whether the 
break-up of heteroclinics on the outside of the jet still occurs. 
By using the Melnikov function from (8), we have been able 
to show that it will occur in a uniform manner all the way down 
to the convergent steering line point. This indicates that the 
region of the convergent steering line presents no more of an 
effective barrier than at higher depths. 

In the numerical simulation, when a second neutral mode 
is added to the streamfunction i.e.,¢ > 0 is set in (9), both the 
heteroclinic and homoclinic tangles are easily observed, as 
seen in Figure 10. The characteristic shape of the heteroclinic 
tangle is easily observable in the figure and it shows that the 
break-up of the heteroclinic orbit results in significant mixing. 
Further numerical experiments also indicate that particles do 
indeed traverse the entire width of the jet. More intricate 
calculations and additional numerical experiments are needed 
to ascertain whether it is a region of more efficacious mixing. 
Of particular significance are simulations for somewhat larger 
amplitudes of the secondary mode, as the dynamics in that case 
is not amenable to perturbation methods like the Melnikov 
method. Related results for two-dimensional maps suggest 
that for such larger perturbations most barriers to transport are 
destroyed and new mixing mechanisms are created. 
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Planforms at the Onset of 
Instability in Double 
Diffusive Convection 
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Introduction 


The subject of double-diffusive convection covers the 
study of fluids in which there are gradients of two or more 
properties with different molecular diffusivities. Many inter- 
esting phenomena have been discovered in the context of the 
ocean, where heat and salt provide the competing properties 
and hence this subject is also known as thermohaline or 
thermosolutal convection. An overview of the area is given in 
the recent articles by Schmitt (1994, 1995). The following 
paragraph takes snippets from them in order to give a quick 
glimpse of the history. 

Over a century before Melvin Stern discovered salt fin- 
gers, W. Stanley Jevons performed the first salt finger experi- 
ment in an attempt to model cirrus clouds. Remarkably, he 
seemed to realize that a more rapid diffusion of heat relative 
to solute played a role in the experiments. However, he over- 
simplified the physics and incorrectly assumed that the “inter- 
filtration of minute, thread-like streams” was a general result 
of superposing heavy fluid over light fluid. Interestingly, Lord 
Rayleigh became aware of these experiments more than two 
decades later. He repeated the Jevons’ experiments at the 
Cavendish Laboratory in Cambridge in April, 1880. The re- 
Sults led him to initiate the study of buoyancy effects in fluids 
by formulating several stability problems for a stratified, but 
non-diffusive fluid. His neglect of diffusion meant that he 


missed an opportunity to discover double-diffusive convec- 
tion. The modern study of double-diffusive convection began 
with Melvin Stern’s article on “The Salt Fountain and Ther- 
mohaline Convection” in 1960. Stommel et al (1956) had 
suggested that a flow (the salt fountain) would be driven in a 
thermally-conducting pipe, but it was Stern who realized that 
the two order of magnitude difference in heat and salt diffusivi- 
ties allowed the ocean to form its own pipes. These later came 
to be known as “salt fingers”. Stern also identified the potential 
for the oscillatory instability when cold, fresh water overlies 
warm, salty water in the 1960 paper. This “diffusive-convec- 
tion” process was demonstrated later by Turner and Stommel 
(1964). 

For the purpose of this article, we will focus on the 
situation with warm salty fluid below cold fluid. The fluid lies 
between two horizontal walls placed a distance L apart, and is 
of infinite extent in the horizontal (x-y) plane. The warm fluid 
tends to rise and the salty fluid, being heavier than fresh fluid, 
tends to fall, so that if we begin with the liquid at rest, then 
increasing the temperature difference will eventually desta- 
bilize the rest state and yield either an oscillatory motion or a 
steady convective motion. In particular, we will address the 
oscillatory onset of instability (the “diffusive-convection” 
problem) and examine a family of patterns that may arise. We 
wish to find new solutions that bifurcate from the rest state, 
and determine their stability to perturbations. 





Methodology 
Stability of Fluid Motions 


We begin by identifying the velocity, pressure, tempera- 
ture and solute concentration that satisfy the governing equa- 
tions and which represent a “base flow”’. The velocity field for 
the base flow that we will investigate in the double diffusion 
problem is zero (the “rest state’’), the temperature and solute 
concentrations depend linearly on the vertical variable z, and 
the pressure field is a quadratic polynomial in z. If this base 
flow is disturbed slightly, then the disturbance may decay, or 
may grow, depending on whether the flow is stable or unstable. 
If the flow is unstable, then the solution will change from the 
base flow to something else, but to what? We will give a partial 
answer by considering small disturbances, and how the distur- 
bances that grow the fastest interact with each other at the 
conditions close to the onset of instability (the “critical situ- 
ation”). We therefore express the total solution as a superpo- 
sition of the base flow and a perturbation, and write down the 
equations that govern the perturbation. In the following, a 
solution will usually mean a solution for the perturbation. 


Spatial Periodicity 


The governing equations for the perturbation to the base 
flow show that if one stands at the origin and looks in any 
direction, the problem looks the same. Many bifurcation prob- 
lems in fluid mechanics involve one or more spatially un- 
bounded directions and a continuum of modes. A full 
description of the set of bifurcating solutions in such a context 
is a rather formidable problem that has not been solved. 
Therefore, one typically imposes some spatial periodicity on 
the problem and then confines attention to solutions satisfying 
this given periodicity. For instance, one may look for solutions 
that are periodic with respect to the vectors 


x ,= W - to ° and x q= W - (0,1,0). 


These vectors are said to span a hexagonal lattice of period W. 
In the x and y directions, the solution is then doubly periodic; 
that is, at any time t, the temperature at the position x = (x,y,z) 
is equal to the value at x+n, x,+n, x, for any integer n,, n,. 
The lattice obtained from this double periodicity is invariant 
under the symmetries of the hexagon;that is, rotation by mul- 
tiples of 60 degrees, reflection across the vectors a; defined by 


4n 4n 13 
#1 Wes 100.02 = WR (— 3-9 Mess Ay 


and reflection across axes perpendicular to the a;. 

The same periodicity condition holds for the other un- 
knowns: the solute concentration, velocity field and pressure 
field. This simplifies the problem immensely and allows us to 
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use a Fourier series expansion of the solution. We next de- 
scribe a method that can be applied to find ordinary differential 
equations for the dominant Fourier amplitudes, and thereby 
find families of spatially-periodic structures on a planar lattice. 
Rolls and hexagons are two of the solutions that are periodic 
on a hexagonal lattice, and often observed in convection 
problems. Other lattices that may be considered are the square 
and rhombus (Golubitsky, Swiftand Knobloch, 1984; Dionne, 
Golubitsky, Silber and Stewart, 1994). 


Center Manifold Theory and 
Dynamical Systems 


The analysis of spatially periodic bifurcating solutions 
and their stability proceeds in two steps. First, the governing 
equations are reduced by projections to a finite number of 
nonlinear ordinary differential equations (Iooss and Joseph, 
1980, Chow and Hale, 1982). These are evolution equations. 
The center manifold theorem is invoked for this purpose. 
Roughly speaking, the theorem says that in the neighborhood 
of criticality, the dynamics is governed by the interactions 
among the finitely many critical modes. This places the prob- 
lem in the framework of dynamical systems. At this stage, 
much of the details of the original partial differential equations 
and boundary conditions now appear in the coefficients of the 
evolution equations. Thus, seemingly different physical prob- 





Figure 1. 


The neutral stability curve for the oscillatory onset (line) and 
steady onset (dashed line) of instability for the stress-free problem 
is shown. The Lewis number ist < 1 
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Figure 2. 


We display the four standing wave solutions: standing rolls, 
standing hexagons, standing regular tnangles, and standing 
patchwork quilt. These are contour plots for an amplitude func- 
tion, e.g., the vertical component of the velocity. (Photographs for 
these figures were taken by Prof. David Wagner, Department of 
Mathematics, University of Houston, Texas.) 
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lems will be seen to reduce to the same types of evolution 
equations. Secondly, the existence and stability of the bifur- 
Cating solutions are determined by analyzing the system of 
ordinary differential equations. This general methodology 
handles a large class of stability problems, where the onset of 
instability is due to just a finite number of modes with the rest 
of the modes being stable, as is typical for viscous flows. 


Governing Equations 


A fluid is assumed to lie between two walls placed at z* 
= 0 and z* = L in the (x*, y*, z*) plane. Asterisks denote 
dimensional variables. The upper plate is kept at a constant 
temperature 6°, with solute concentration S”, and the lower 
boundary is kept at a higher temperature O*, and a higher 
solute concentration S*,; 6°, < @*,,S*, <S, - At temperature 
9°,, the fluid has thermal coefficient of volume expansion a, 
solute coefficient of volume expansion 6, thermal diffusivity 
Ky, Solute diffusivity K,, viscosity |1, density Pp, and kinematic 
ViSCOSity V=[1/Pp. 

Dimensionless variables (without asterisks) are as fol- 
lows: (x, y, z) = (x*, y*, z*) /L, t= Kpt*/L?, v = v*L/k;, 0 = 
(8*-8*,)/(8*,-0*,),S= (S*-S*,)(S*,-S*,), Pp = p*L7/(Poky). 


A thermal Rayleigh number R,, a salinity Rayleigh number 
Rg, a Prandu number P, and Lewis number t are defined by: 
*,- 0 )L’/(x;y), 


R,= ga (6 R,=gB(S | -S ) L/(x,»), 


P= V/Kr, t= Ks/ Kr. 
where g denotes gravitational acceleration, and a dimension- 
less measure of gravity is G=gL3/k,’. 

A basic solution is 


v=0, 6=1-z, S= 


p=py- Gr + (Ry- Ry) P(e), 


1-—z, 


where py is a constant. Since we are concerned with the 
stability of this rest state and with solutions bifurcating from 
it, we set up the problem for the perturbations, 5, p and v = 
(u, v, w) to this solution. The equations satisfied by these 
variables are: 


the transport equations for the temperature and solute concen- 
tration 


§ -w-Ab=-(v-V)G, (2) 


3 —-w-taS=-(v-V)S, (3) 


the Navier-Stokes equation with the Oberbeck-Boussinesq 
approximation 


v—PAv+Vp+(RsPS —R7P8) e,=-(v-V)v, (4) 


and incompressibility 
div v=0. (5) 


Values for the temperature and solute concentration are pre- 
scribed at the walls: 


G=0, 3=0, a 2z=0 and z=1. 


The velocity boundary condition at the walls is the “no-slip” 
condition: 


v=0, at z=0 and z=1. (7) 


Alternatively, one may use the “stress-free” boundary condi- 
tion (zero normal velocity and zero shear stress): 


, ow 


> (8) 


For the stress-free problem, the eigenfunctions are available 
in closed form in terms of trigonometric functions (Baines and 
Gill, 1969, Huppert and Moore, 1976, Tumer, 1979, Nagata 
and Thomas, 1984, Knobloch, 1985, Deane et al, 1987, 
Knobloch et al, 1986), and therefore the weakly nonlinear 
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Figure 3. 


We display the three travelling wave solutions. The upper pattern 
shows travelling rolls. The lower patterns are the travelling patch- 
work quilt (1) and the travelling patchwork quilt (2). The arrows 
indicate the direction of travel. (photographs by David Wagner) 
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analysis can be done explicitly. This explicit representation is 
not available for the case of no slip, and a computational 
approach is then used for the analysis. The answer to the question 
of which type of velocity boundary condition is pertinent, lies in 
which phenomenon one is trying to describe. The use of the 
stress-free condition is relevant for modeling situations where the 
fluid layer is not bounded by walls but by free surfaces; for 
example, double diffusive actions may be taking place in a 
sandwiched layer at some depth say, in the ocean. 

The set of unknowns, v, p,5 and § and §, is denoted by ®. 
Equations (2) - (8) can be written in the schematic form 


L® =N, (®,®), 


Here the operator L incorporates the linear terms, while N, 
stands for the quadratic terms. 


Linear Stability Analysis 


For the case of infinitesimally small perturbations, we 
linearize (9). The rotational symmetry of the problem about 
the vertical axis shows that we can ignore the y dependence 
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and pose the linearized problem in the x-z plane. The method 
of separation of variables is used to express the solutions as 
exp(iax+o) multiplied by a function of z. The real linear 
operator L has the form A +4, so that we are seeking eigen- 
functions € which satisfy (A+ o B)C = 0. The wavenumber of 
the perturbation is & and the growth rate is Re o. We wish to 
map out the regions of stability (Re o < 0) and instability, by 
locating the neutral stability curve (Re o = 0). 

Figure 1 shows the neutral stability curve for the stress- 
free problem (Nield, 1967, Baines and Gill, 1969) for the case 
t < 1, which is typical for many applications, including the 
ocean media. The neutral stability curve consists of two lines 
meeting at a Takens-Bogdanov point. The steady onset 
(dashed line) is along 


Rs 2744 
os tel a 


and the oscillatory onset is along 


P+t t, 27x‘ 
rr =(F isa +1) (1+5) 4 


The critical wavenumber is a= 1/ V2, the same along the entire 
neutral stability curve. 

For the no slip problem, the onset conditions need to be 
computed numerically. For our numerical results, we use the 
Chebyshev-tau method, which approximates the eigenvalues 
o with infinite-order accuracy. This is a spectral method 
(Orszag, 1971, Gottlieb and Orszag, 1983) in which the z-de- 
pendence of the eigenfunctions is expressed in terms of Che- 
byshev polynomials. In this manner, the linear stability 
problem is discretized into a matrix eigenvalue problem. The 
neutral stability curve for the steady onset is 


a = 3.12. 


Rs 
Rr= < + 1707.765, 


The neutral stability curve for the oscillatory onset is more 
complicated. This curve lies above that of figure 1, in that the 
critical thermal Rayleigh number is higher at each R,. More- 
over, the critical wavenumber increases with Rg. 


Nonlinear Analysis 
Background 


In the two-dimensional case, there is a critical eigenfunc- 
tion representing waves traveling to the right (with wavenum- 
ber a) and there is another one traveling to the left (with 
wavenumber -c). These interact nonlinearly to produce trav- 
eling and standing wave solutions. In this situation, if both the 
standing and traveling waves bifurcate supercritically, then 
one of them is stable (Ruelle, 1973). Past work on nonlinear 








Figure 4. 


We display the oscillating triangles. Upper patterns: t=0, t=T/24, 
T is the period. Middle patterns: t=T/12, t=T/8. Lower pattern: 
t=T/6. (photographs by David Wagner) 
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analysis has been on the stress-free problem (Knobloch, 1985). 
Deane et al (1987) derive the amplitude equations for the left- 
and right-traveling waves. At the least, they need to calculate 
up to third order in the amplitude functions. The calculation of 
the coefficients in the equation yields the Landau coefficient 
which determines whether the solution is supercritical or not. 
They find that the third-order Landau coefficient for the trav- 
eling wave turns out to be zero. This degeneracy appears to be 
a coincidental one (Knobloch et al, 1986) due to the adoption 
of the stress-free condition. The standing wave solution is 
supercritical, indeterminate and subcritical, respectively, for 
Rg less than, equal to and greater than 14585.4 for the case 
t=10'* P = 1. At the fifth order, Deane et al find that the 
traveling wave solution bifurcates supercritically, so that there 
is a range of R, for which the traveling wave solution is stable. 


What are the third order Landau coefficients, and what is 
the stable solution for the no slip case? Computations at the 
above values of Lewis and Prandtl! numbers for R.=10* 
14585.4, and 10° reveal that the solutions are both unstable 
(Renardy, 1993). As expected, the selection mechanism is 
sensitive to the boundary conditions. 

What other periodic solutions are there in the fully three- 
dimensional planform problem, and are they stable to the more 
general doubly periodic disturbances in the x-y plane? The 
two-dimensional solutions we have discussed correspond to 
standing rolls and traveling rolls in three dimensions. The 
three-dimensional stress-free problem has been studied by 
Nagata and Thomas (1984). For oscillatory onset, they present 
the linear stability of standing rolls, standing squares and 
standing hexagons to perturbations restricted to have the same 
symmetry as the solutions. For example, they examine 
whether hexagon pattern solutions are locally asymptotically 
stable with respect to hexagon pattern disturbances, and 
whether the roll solutions are stable with respect to roll-like 
disturbances. The next issue is the competition, say between 
rolls and hexagons, and the investigation of the traveling wave 
solutions. 


Center Manifold Theorem and 


Dynamical Systems 


The fluid properties are regarded as fixed and we define 
a bifurcation parameter A=R;-Ryc, where R, is the thermal 
Rayleigh number and Ry, is its critical value. Solutions ® with 
the periodicity of the hexagonal lattice are sought. Details of 
this calculation are found in Renardy (1993). 

At criticality (A=0), there is a pair of complex conjugate 
eigenvalues + iw of the linearized problem, and each of them 
has six eigenfunctions. For A near 0, -41(A) denotes the eigen- 
value which arises from perturbing -iw. We denote by C,(A), 
k=1,2,...,6 the eigenfunctions belonging to -1(A), i.e. 


L(- w(A))C (A) =0, 


and those belonging to - ji() are their complex conjugates. We 
will need to compute the eigenfunctions at A=0. 

The critical eigenfunction computed in the (x,z)-plane has 
the form C, =¢(z) exp (ia, - x), where a,=(a,0,0), x =(x,y,z), 
and « is the critical wavenumber. We may visualize the six 
eigenfunctions as follows. The vectors + a,, + a, and + a, 
defined by equation (1) emanate from the center of a hexagon 
and terminate at its six vertices. The critical eigenfunctions are 
waves propagating in the directions of a,, a,, a, and - a), - a, 
- a,. We denote the (complex) amplitude of the wave propa- 
gating in the direction of a; by z, and the amplitude of the wave 
propagating in the direction of - a; by z;,,. 

We can decompose a solution ® in the form 


O=0, +, 
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6 


k=1 


where z, are the complex-valued amplitude functions and ‘¥ 
represents a superposition of eigenvectors (and possibly gen- 
eralized eigenvectors) belonging to stable eigenvalues.To use 
the Center Manifold Theorem, we need to have the fact that ‘VY 
consists of stable modes. The expression on the right hand side 
of (10) converges to ® in the mean square. The ¢ and the 
functions in ‘¥ form a complete set. We wish to find equations 
that govern the evolution of the amplitude functions z,. In order 
to do that, we need to find an explicit approximation to ‘Y for 
small perturbations; at this stage, we invoke the Center Mani- 
fold Theorem. 

The Center Manifold Theorem states that, in a neighbor- 
hood of ®=0, there is a manifold (called the center manifold) 
of the form ‘¥=1(Z, ,2>,23,24,Zs.Z,) with the following proper- 
ties: 


l. All solutions with initial data on the center manifold 
remain on the center manifold as long as they remain 
small. 


All small periodic solutions lie on the center manifold. 


The stability of a small periodic soluuion is determined 
by its stability within the center manifold, in other 
words, all Floquet multipliers corresponding to direc- 
tions outside the center manifold are stable. 
These properties imply that only quadratic terms in the asymp- 
totic approximation to the center manifold are required. Thus, 
the first property above is used to write the following expres- 
sion 


6 


Y,=2 Re ( bx 2; Zz Wij + 2; 2; Xi)» 
iy=1 


Y=, +..., 


where the dots indicate terms of higher than quadratic order. 
The y;; and x;; contribute to second harmonics and the mean 
flow mode. 

The second property above says that the solutions of 
interest (the small periodic solutions, taking into account the 
nonlinear terms) are found on the center manifold. These are 
solutions close to the linear eigenfunctions. There need not be 
a center manifold if amplitudes become large. 

By forming appropriate inner products and projections 
onto appropriate spaces, we are able to calculate the functions 
yj; and x;;. The total solution ®, accurate to second order, is 
now substituted back into the governing equations, leaving the 
z, as the only unknowns. After further projections, using the 
symmetries of the hexagonal lattice and the theory of normal 
forms, the following reduced system is obtained: 
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Figure 5. 


Wavy rolls (1). Top first row: t=0, t=T/16, T is a period. Second 
row patterns: t=T/8, t=3T/16. Third row: t=T/4, t=5T/16. Fourth 
row: t=3T/8, t=7T/16. Fifth row: t=T/2. (photographs by David 
Wagner) 
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Figure 6. 


Tusted patchwork quilt. Upper patterns: t=0, t=1/24, Tis a period. 
Middle patterns: t=T/12, t=T/8. Lower pattern: t=T/6. (photo- 
graphs by David Wagner} 





t=V/12 





2 ; 
Ge + Fi @122,23,2402526A) = 0,i = 1,....6, (11) 


Fy (21 29,25 :Zqs25.2gA) =p (A) 2, +O, (A) Iz, Pz, +0, (A) 12,7 2, 
+04, (A) (1 zy? +125 I) 2, + 01g (A) 
(125 PF +126 I) 21 + 5 (A) (29 25 + 23 26) Hy t-- 
The dots denote terms of higher order, and c,,...,0, denote 
Landau coefficients. The F,, i=2,..,6 are related to F, by prop- 
erties of the hexagonal symmetry. We have now reduced the 


governing equations to a system of six complex nonlinear 
ordinary differential equations. 


Pattern Selection on the Hexagonal 
Lattice 


In (11), we set u = (2, ,2,,23,24,25,Z,), and obtain 


& + F (wd)=0, (12) 
where uw eIR' and A is a real parameter. Starting with an 
equation of this general form, Roberts, Swift and Wagner 
(1986) use group theory to classify new solutions and deter- 
mine their stability properties by considering their symme- 
tries. This approach has allowed them to find a larger family 
of solutions than was previously available. Their results have 
a generality in that they apply to other problems that can be 
expressed in the form of (12) and have the same symmetries. 
They find eleven qualitatively different classes of bifurcated 
solutions: (i) Standing rolls, (ii) Standing hexagons, (iii) 
Standing regular triangles, (iv) Standing patchwork quilt, (v) 
Travelling rolls, (vi) Travelling patchwork quilt (1), which is 
always unstable, (vii) Travelling patchwork quilt (2), (viii) 
Oscillating triangles, (ix) Wavy rolls (1), (x) Twisted patch- 
work quilt, (xi) Wavy rolls (2). These are solutions which have 
enough symmetries so that the twelve degrees of freedom are 
reduced to two. For example, invariance under rotations by 
multiples of 60 degrees and reflections across the axes of the 
hexagon yields a two-dimensional solution set where the zi are 
equal, called the standing hexagon. The conditions for stability 
of each solution with respect to perturbations with the double 
periodicity is also provided in Roberts et al. To apply these 
results to a specific problem, it is necessary to calculate the 
coefficients that appear in (11). 

Figures 2 - 7 are contour plots of the vertical component 
of the velocity for each of the solutions. Solutions (i) - (iv) are 
standing waves, i.e., the patterns simply oscillate with time, 
and are shown in figure 2. These plots can be thought of as 
pictures taken at time zero. Solutions (v) - (vii) are traveling 
waves, i.e., the pattern moves to the right without changing its 
shape, and are shown in figure 3. The four remaining solutions, 
shown in figures 4 - 7, are neither standing nor traveling waves. 
The time dependence is indicated by a series of pictures. For 
example, figure 4 shows the oscillating triangle at t=0, 1/24, 
1/12, 1/8 and 1/6 of a period. The pattern evolves through the 
rest of the period before going back to the t=0 picture. 


Numerical Results for the Oscillatory 
Onset for Double Diffusion 


For Lewis number t=10"'” and Prandt! number P = 1, we 
have seen that traveling waves are stable in the two-dimen- 
sional stress-free problem for a range of solutal Rayleigh 
numbers Rg, while it is unstable for the no slip problem. We 
pursue the no slip problem as a model for the fluid bounded 
by walls, and examine the stability of the 11 solutions in the 
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three-dimensional case. To achieve an oscillatory onset, the 
solutal Rayleigh number should be larger than approximately 
10°; results in the range 10°< R,<10® will be presented. It is 
found that the solutions are unstable until the solutal Rayleigh 
number becomes relatively large. As the R, increases, so does 
the R, at criticality, since the solute concentration effect must 
counteract the thermal effect to provide for the oscillatory 
onset (at Ro=2 x 10°, the critical R, is of order 1 million). For 
solutal Rayleigh numbers between 8 x 10° and 8.5 x 10°, the 
standing rolls are stable and the other 10 solutions are unstable. 
There is a window of R, from 9.0 x 10° to about 2.0 x 10°, 
where all 11 solutions are unstable. For R,~ 4.0 x10°, either 
the standing hexagons or the standing regular triangles is 
stable but not both, and the remaining 9 solutions are unstable. 
As a reminder of the scales involved, the typical properties for 
water (Drazin and Reid, 1981) are: & = 5 x 10*K"', 0,” - 6,” 
< 10 K, x, = 10° cm/sec, v = 10 cm/sec, g= 10° cm/sec’. 
Substituting these into the definition for the thermal Rayleigh 
number to be 1 million, we find L’= 2, L being measured in 
cm. 

The heating of cold salty water is more appropriately 
modelled by a Prandtl number of P=10.0 and Lewis number 

=0.01. Numerical results were obtained for solutal Rayleigh 
numbers up to order a million. At low Rayleigh numbers, the 
solutions were found to be unstable. For R, between 10° to 
10’, the standing rolls are stable and the other solutions are 
unstable. For R,~ 2x 10’, three of the standing wave solutions 
are found to be stable. These are the standing rolls, the standing 
patchwork quilt, and either the standing hexagons or the 
Standing regular triangles. If this situation were realized in 
nature, we may see a coexistence of these solutions, a patch of 
one lying next to another. 

We have obtained numerical results for two sets of Prandtl 
and Lewis numbers; the 11 solutions tend to be unstable for 
low Rayleigh numbers, while a number of standing wave 
solutions may be stable at higher Rayleigh numbers. 


Outlook for the Future 


The center manifold reduction method is applicable to 
onsets of instability, whether steady or oscillatory or both, in 
many different viscous flow situations, and continues to be 
applied to a variety of problems. The group theoretic approach 
for finding classes of solutions to the nonlinear ordinary 
differential equations by taking advantage of the symmetries 
in the problem is particularly useful for convection problems; 
recent examples include ferrofluids (Silber and Knobloch, 
1991), viscoelastic liquids (Renardy and Renardy, 1992), in- 
terfacial stability in a two-layer system (Renardy and Renardy, 
1988; Joseph and Renardy, 1993) and magnetoconvection 
(Clune and Knobloch, 1994). On the one hand, these analytical 
results serve to explain observed phenomena, and on the other, 
they encourage future search for these patterns in the labora- 
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tory. For instance, we have found that three different patterns 
can coexist in certain parameter ranges. It would be interesting 
to see which one(s) of these can be set up under experimental 
conditions. 

This article has taken the oscillatory onset for the double 
diffusion problem and highlighted the dynamical systems 
approach coupled with group theory. Apart from the oscilla- 
tory onset, the double diffusion problem can have a steady 
onset or the case where the oscillatory onset merges into a 
steady onset. In particular, the case of the steady onset is 
relevant for the investigation of salt fingers. The subject of salt 
fingers, and their stability, is an ongoing area of research, 
where mathematical modeling together with the type of analy- 
sis presented here provide fertile grounds for future interdis- 
ciplinary collaboration. For instance, why are square patterns 
commonly observed (Proctor and Holyer, 1986)? and why are 





Figure 7. 


Wavy rolls (2). Upper patterns: t=0, t=7T/24, T is a period. 
Middle patterns: t=1/12, t=T/8. Lower pattern: t=T/6. (photo- 
graphs by David Wagner) 











there triangular and asymmetric patterns in certain regions of 
the ocean (Schmitt, 1994b)? 
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